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Abstract 

We study numerically the Coulomb interacting two-particle stationary states of the Schrodinger equation, where 
the particles are confined in a two-dimensional infinite square well. Inside the domain the particles are subjected 
to a steeply increasing isotropic harmonic potential, resembling that in a nucleus. For these circumstances we have 
developed a fully discretized finite difference method of the Numerov-type that approximates the four-dimensional 
Laplace operator, and thus the whole Schrodinger equation, with a local truncation error of 0{h 6 ), with h being 
the uniform step size. The method is built on a 89-point central difference scheme in the four-dimensional grid. As 
expected from the general theorem by Keller [Num. Math. 7, 412 (1965)], the error of eigenvalues so obtained are 
found to be the same order of magnitude which we have proved analytically as well. Based on this difference scheme 
we have obtained a generalized matrix Schrodinger equation by vectorization. In the course of its numerical solution 
group theoretical methods were applied extensively to classify energy eigenvalues and associated two-particle wave 
functions. This classification scheme inherently accounts for the symmetry property of the two-particle state under 
permutation, thereby making it very easy to fully explore the completely symmetric and antisymmetric subspaces 
of the full Hilbert space. We have obtained the invariance group of the interacting Hamiltonian and determined 
its irreducible representations. With the help of these and Wigner-Eckart theorem, we derived an equivalent block 
diagonal form of the eigenvalue equation. In the low-energy subspace of the full Hilbert space we numerically computed 
the ground state and many (m 200) excited states for the noninteracting as well as the interacting cases. Comparison 
with the noninteracting exact results, which we have found analytically too, reveals indeed that already at a modest 
resolution of h m 1/15 our numerical data for the eigenvalues are accurate globally up to three or more digits of 
precision. Having obtained the energy eigenvalues and eigenstates we calculated some relevant physical quantities 
with and without interaction. These are the static two-particle densities labeled by irreducible representations, one- 
and two-particle density of states and some different measures of entanglement, including the reduced density matrix 
and the von Neumann entropy. All these quantities signal concordantly that the symmetric states (with respect 
to permutation of particles) as well as their energies are affected by the interaction more dramatically then their 
antisymmetric counterparts. 
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1. Introduction 

The importance of numerical solution of the Schrodinger equation for analytically unsolvable potentials 
can hardly be overestimated in any fields of modern physics. In one-dimension (ID) the search for either 
exact or numerical solutions with ever increasing accuracy has been in the forefront of research for many 
decades now 1 1 121 . One of the reasons for the ID case has earned that much attention lies partially in the fact 
that many higher dimensional problems have spherically symmetric potentials which then allows separation 
of variables and the study of the radial equation only. Numerical methods for obtaining bound state as well 
as scattering state solutions can be broadly classified into two types: numerical integration or in other words 
the shooting methods [314] and matrix methods [2[5 6 7 8 9 10 . Also, there are analytical approaches too 
like perturbative treatments or the Rayleigh-Ritz variational principle, see Ref. [6J and references therein. 

The fact that this fundamental equation of non-relativistic quantum physics is a linear second order 
differential equation without a first order term makes it especially suitable to be studied by, beyond the 
classical finite difference methods |5I6) , Numerov's method as well |6lllll2ll3ll4j . As to the classical methods, 
these are the most simple finite difference schemes based on the fact that the differential operator representing 
the kinetic energy can be expressed as a series of central difference operators [5 15 . Truncating this expansion 
at the first, second or third terms yields algebraic eigenvalue equations with tri-, penta- or heptadiagonal 
matrices, respectively. Also, it turns out that the local truncation error characterizing the approximation 
of the continuous differential equation with the difference equation is, in these cases, 0(h n ), with n = 2, 
4 and 6 and h is the uniform step size. More interestingly, it can be proven analytically that the errors of 
eigenvalues so obtained have the same order of magnitude as the truncation error [16) . Practically it means 
that in the heptadiagonal case for example seven to eight digits of precision can be easily reached. Because of 
this high accuracy they have been used extensively to calculate the properties of the ground state and some 
low energy excited states of many different potentials. These include the anharmonic oscillator potential, 
the symmetric double well, the Razavy potential and many others [BJ. 

Despite the success of the classical finite difference schemes mentioned above, the most widely used method 
for solving the Schrodinger equation is, however, the Numerov method [6|ll|12|13|14j . This is because the 
first approximation already yields very accurate eigenvalues with error at most 0(h A ). This is also true in two 
and higher dimensions, though the procedure is slightly complicated by the fact that the recurrence relation 
involves second-neighbors as well and so some additional information or physically motivated approximation 
is needed close to the boundary in order to start the algorithm [91 14] . Though its performance, in the 
sense of accuracy, is obviously not as good as the classical heptadiagonal approximation, with higher order 
approximations for the derivatives |13|17 18 or with extrapolation techniques [7 8 9J (Richardson, Aitken) 
it can be refined easily up to 0(h 6 ) or higher. 

As to higher dimensions, it is rather interesting that despite the availability of a great deal of numerical 
methods, developed primarily for ID, only a few of them have been adopted to two or higher dimensional 
Schrodinger problems. Perhaps the most natural and simple approach of all is the classical fully discretized 
five-point method with an error of 0(h 2 ) [7|8|19j . This was successfully applied in two-dimensions (2D) 
for potentials that lack rotational symmetry. Later both the standard Numerov method and its extended 
version have been transported to 2D in order to obtain, respectively, 0{h A ) and 0(h 6 ) accurate results 
for the lowest eigenvalues of some Coulombic potentials [S]. Also, partial discretization methods have been 
proposed recently for the study of the highly anisotropic Henon-Heiles potential [10] . 

Our aim in this paper is to study numerically an interacting two-particle problem in two-dimensions. The 
geometry of the problem and the number of particles imply naturally that we have to face a four-dimensional 
Schrodinger equation. The method we apply is the full discretization of the differential equation which we 
solve with a generalized Numerov-type approach. Generalized in the sense that (i) the algorithm is 0(h 6 ) 
accurate, and (ii) to our knowledge Numerov's method has not been ported to dimensions greater than two 
[20j . From the physics point of view we are of course primarily interested in the study of interaction and so the 
mathematical background we have developed along the way could be considered as a necessity. Nevertheless, 
we feel that it is alone very interesting because it reveals the connection of the high-dimensional finite 
difference method with the theory of matrices composed of commuting blocks. It could be utilized elsewhere 



2 



and in other dimensions as well. The paper is organized as follows. In Section [3 we start with the precise 
formulation of the quantum mechanical problem, then develop the necessary algebra. In Section [3] we study 
the group theoretical aspects of the problem. Sections [4] and [5] are devoted to our numerical results without 
and with interaction, respectively. Finally, Section [6] summarizes our conclusions. 

2. Schrodinger equation and stationary states 

2.1. Problem formulation 

Let us start by the formulation of the quantum mechanical problem of two identical interacting particles. 
The quantum theory we are about to explore is non-relativistic and therefore the spin degrees of freedom 
of particles will not show up explicitly in the Schrodinger equation. For the same reason neither spin-orbit 
coupling will be considered here. The total two-particle wave function is thus separable as ^(ri, o~x; t 2 , cr 2 ) = 
"0(ri, r2)x{< J ii where \ is the spin dependent component for the construction of which the knowledge of 
statistics (Fermi or Bose) and the appropriate spin quantum number is essential, while ip is the orbital part 
governed by the usual time independent Schrodinger equation 

(-|^(Ai + A 2 ) + C/(n) + U(r 2 ) + y(|n - r 2 |)^ V(ri,r 2 ) = £V(ri,r 2 ). (la) 

Here U is the external potential that acts on both particles separately and V is the pair interaction de- 
pending only on the distance between them. In order to determine possible energy eigenvalues E and the 
corresponding states ^ unambiguously, that is the spectrum of the Hamiltonian, one needs to supplement 
Eq. (jlap with a boundary condition. In the paper we consider a two-dimensional square shaped potential 
well centered on the origin 

D 2 = [-b,b]x[-b,b], (lb) 
inside of which U(r) is a well-behaved smooth function, but outside U(r) = Uq — > 00. This infinite potential 
barrier gives a constraint on the wave function ipiji, r 2 ): for both arguments on the boundary 

1>\bd, = 0. (lc) 

Now that the dimensionality of the problem is specified, the Laplacian in Eq. (fTa|) is expressed as = 
d 2 /dx 2 + d 2 /dy 2 , i = 1,2. Clearly, it is the Descartes coordinate system that best suits the rectangular 
domain. From purely mathematical point of view Eqs. ([I]) define a (self-consistent) boundary value problem, 
where the parameter E has to be determined as well. It is its precise numerical solution that is the primary 
object of the present investigation. 

2.2. Method of errors 0{h 6 ) 

From now on, to facilitate numerical calculations, we use atomic units where the unit of length is the 
Bohr-radius ao — h 2 /(me 2 ) and energy is measured in the Rydberg unit i?R y d = e 2 /(2ao) = h 2 /(2mafy. It 
is also convenient to introduce unified notation for space variables as 

(21, 1/1, £2, y 2 ) -> (x 1 ,x 2 ,x 3 ,x 4: ) = x, (2) 

and an effective potential with 

U(x) = U(r 1 ) + U(r 2 ) + V(\r 1 -T 2 \). (3) 

Equation ^ shows that in the unified four-dimensional configuration space the first two dimensions com- 
posed of (xi,x 2 ) belong to one particle, say the "first", and the complementary subspace to the other, the 
"second" . With all these the Schrodinger equation can be reformulated as a Dirichlet problem for Poisson's 
equation 

AV-(x) = -/(x), (4a) 
/(x) = (E - ^(x))V(x), (4b) 
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where A = J^.- L d 2 /dxf, the region of confinement is a four-dimensional (4D) hypercube 

D 4 = [-b, b] x [~b, b] x [-6, 6] x [-6, 6], (4c) 

and the boundary condition in Eq. (jlcj) now reads as 

*P\dD 4 = 0. (4d) 

It is clear that Eqs. ((4]) can also be thought of as a one-particle Schrodinger equation for an abstract four- 
dimensional particle that is confined in a box where the "external" potential is given by Eq. ^ . This picture 
will turn out to be very useful when we turn our attention to the representation theory of the invariance 
group of the Hamiltonian in Section [3] Before that, however, we develop the necessary algebra: based on a 
generalized Numerov-type approach we derive an equivalent matrix Schrodinger equation that provides the 
wave function on a finite grid [619110). To this end we divide the interval [—6, b] into (n + 1) equal segments 
with increment h given by h — 2b/ (n + 1). By this procedure a linearly spaced 4D cubic grid is obtained, 
where the pth lattice point in the ith dimension is Xi iV = —b + ph, p = 1, . . . , n and i = 1, . . . ,4. Further, we 
use the simplified notation 

v{xi tP ,X2 t i,X3,k,X4,i) = VpiM, p,i, k,l= 1,2, ...,n, (5) 

for any multivariable function v. Now define first-, second-, third- and fourth-neighbor difference operators, 
respectively. These are the natural generalizations of the central difference operators (three-point, five- 
point, . . . ) applied commonly in one-dimension |5|6j and two-dimensions |7|8|9|19j . In 4D the first-neighbor 
difference operator reads 

Oi4>pikl = ^ {lpp+a,ikl + 1pp,i+a,kl + i ] pi,k+a.l + 1ppik,l+a) ~ S^pikl- (6a) 
a=±l 

On the right hand side g\ = 8 is the total number of first-neighbor sites. For clarity we shall note here that 
□ i acts on the function and not on its values, thus on the left hand side {\3iip} p iki would be more precise. 
Nevertheless, we stick with the shorter notation as there should be no confusion by dropping those extra 
braces. Similarly the second-neighbor difference operator is 

02lppikl — ^ {4'p+a,i+l3.kl J ri'p+a,i,k+li.l +4>p+a,ik,l+f3 +4>p,i+a,k+P,l +lpp,i+a,k,l+P +lppi,k+a,l+p) — 24-0 p jfc; , 
a,P=±l 

(6b) 

where = 24 is the total number of second-neighbor sites in 4D. In the third-neighbor case 

Oslppikl = (VV+a,i+/3 : fc+7,i + 1pp+a,i+P,k,l+~/ + VV+a,i,fc+/3,i+7 + i>p,i+a,k+P,l+l) ™ 32?/>pifci, (6c) 

a,/3,7=±l 

with g^ — 32, and finally 

t^ilppikl = 1pp+a,i+P,k-H,l+5 ~ 16l/Vifci + [Oi~lp p ikl}h^2h- (6d) 

a,/3,7,<5=±l 

Note that the fourth-neighbor sites of a given lattice point are at a distance \J\h — \J a 2 + /3 2 + 7 2 + S 2 h, 
so they can be situated along the 4D space diagonal as well as along the coordinate axes. Hence, there are 
altogether g 4 = 24 of them. This fact is represented by the last term in Eq. (|6d[) . With these expressions 
at hand we can make a linear combination with arbitrary constants a, [3, 7 and 5, and a subsequent Taylor 
expansion in h up to order h 6 yields 

(oOi + /3D 2 + jn 3 +6n 4 )ip pik i = h 2 {a + 6/3 + 12 7 + 125)Ai/> 

+ ^([« + 6/3 + 12 7 + 245] + • • • + ^4) + [12/3 + 48 7 + 485]^^ +...)) 
+ jj^([a + 6/3 + 12 7 + 725] (V x? + • • • + + [360 7 + 7205] (V x?x | xi + . . .) 

+ [30/3 + 1207 + 1205](^ X 2 X 4 + i> <xl + ^ ?x | + ^ <x% +...)). (7) 
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Here, all derivatives on the right hand side must be evaluated at the lattice point p, i, k, I. As to the partial 
derivatives, we used shorthand notation, for example 



{mm'")*' (8) 



The first term on the right hand side of Eq. (J7JI is of order h 2 and consists of only second derivatives. 
Similarly, the second (third) term is of order /i 4 (/i 6 ), consisting of only fourth (sixth) order partial derivatives. 
Moreover, the orders of the derivatives with respect to each variable are even. The fact that this must be so 
and that there are no terms proportional to odd powers of h are justified by the definitions of Namely, 
all the neighborhoods of a given lattice point appearing in Eqs. ^ are in a sense complete: they are all 
symmetric under 4D spatial inversion and reflections, thus odd terms must drop out indeed. With simple 
combinatorial arguments it is easy to see that there are 6, 4 and 12 terms exhibiting structure like ^ X 2 X 2, 
%j) x 2 x 2 x 2 and ip x 2 x i, respectively. Equation ([7]) holds for every analytic function and so are 

A 2 ip = V> x 4 H h ip x * + 2(i; x 2 x 2 + ip x 2 x 2 + ...), (9a) 

A 3 "0 = 1p x 6 H h 1p x 6 + 3{i' x 2 x 4 + 1p x 4 x 2 + 1p x 2 iX i + l{j x i x 2 +...)+ 6(^2^2 +...), (9b) 



and 



(d x 2 x 2 + • • • + 8 X 2 X 2 )Alp = 3(lp x 2 x 2 x 2 +...)+ Ipxjx^ + ^ x\x\ + ^ x\x\ + ^ x\x\ + • ■ ■ ■ (9c) 

These identities can be verified easily by direct calculation. 

So far the calculations might seem rather artificial and the formulas apply to all differentiable function 
including the true (yet unknown) physical solution of Eqs. @. Now, the crucial recognition that validates 
prior work is the fact that one can choose the parameters of Eq. (J7JI in such a way that its right hand side 
is a differential form of Aip only. It is easy to see that the term of order h A is much like Eq. (|9"a|) , one only 
has to choose the prefactors appropriately. In a similar manner one can also realize that the term of order 
h 6 consists of derivatives that appear in Eqs. (|9b|) and (|9c|). Therefore, with a suitable choice for a, (3, 7 and 
5 we (possibly) can achieve that to be proportional to 

A 3 ^ + v(d xjx 2 2 +--- + d x 2 sX 2 4 )AV, (10) 

where rj is introduced as a fifth unknown. With all these findings a system of linear equations is obtained 

2 (a + 6/3 + 12 7 + 24(5) = 12/3 + 48 7 + 485, (11a) 

3(a + 6(3+ 12 7 + 72(5) = 30(3 + 120 7 + 1205 - 77, (lib) 

6(a + 6/3 + 12 7 + 72(5) = 360 7 + 720(5 - 3rj, (11c) 

with solution 

a = 12 7 , (12a) 

13 = 7 + 8(5, (12b) 

r) = 6O7. (12c) 

As there are only three equations it is not surprising that among the five unknowns two, say 7 and 6, 
remained undetermined. The question how to fix these degrees of freedom is very interesting from both 
purely mathematical and physical points of view. We will consider this issue in detail in Subsection 12 . 31 and 
we will see that it is rather nontrivial and requires very careful analysis. Now, putting all these expressions 
together and remembering the fact that the physical solution satisfies Eq. ()4a|) we end up with 

- ^2 ( 12 7 n i + (7 + 8<5)D 2 + 7^3 + Sn 4 ) tppiki =w p m, p,i,k,l = l,...,n, (13a) 

where 

w pM = (7 + 2S)f pikl + h 2 + (A f) inkl + ((^0 + if) ( A2 ^ fci + + ■■■ W) ■ ( 13b ) 

This is the equivalent of Poisson's equation, now discretized on a 4D grid [19j . If we were to solve a true 
Poisson's equation numerically with a given differentiable function /, this would be the starting point. Then 
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of course we could keep the formalism as easy as possible by taking S = and 7 = 1, without affecting 
precision. However, in our case / is unknown as it is itself determined by if>, see Eq. (|4b() . It follows at once 
that this equation in the present form is not (completely) a linear difference equation, because w P iki still 
involves continuous derivatives. This difficulty, however, can be easily overcome with the technique developed 
so far: one has to express the right hand side of Eq. (|13bj) as a difference operator acting on /, with the 
constraint that the error so introduced must not exceed 0(h 6 ). This is because the local truncation error 
(we have already made by cutting the series expansion in Eq. ([7])) on the right hand side of Eq. Q13ap is at 
most 0(h 6 ). Making use of Eq. (0) again, now written for a function / and the requirements regarding the 
error of discretization, a second system of linear equations is found 

a' + 6/3' + 12 7 ' + 126' = J- + -, (14a) 

12 5 

a' + 6/3' + 12 7 ' + 246'= ^ + ^, (14b) 



12/3' + 48 7 ' + 485' = 2 ( ^ + +— , (14c) 



whose solution is 



30 15' 

. 30 I / ! " 

a' = 127'-^, (15a) 
* = + * (15b) 

5 ' = -2i0-l4- (15C) 
As in Eqs. (|11[) before, the number of equations is again less than that of unknowns (with primes), and as 
such one variable, say 7', varies freely. Putting these results together leads finally to 

1 



(127Cll + ( 7 + 8<5)D 2 + 7CI3 + SO*) iflpikl 

26 „ A m , ,rn ( 7 , S 



30h 2 

= (, + 2S + (127' -£) Dl+ ^ + »- 47') + 7'Ds - + ^ ) D 4 J Uu (16) 

where, as required, the error is at most 0(h 6 ). This is again a discretized version of Poisson's equation that 
provides ippiki given the set of fpiki is known. Though its solution is the same as that of Eqs. (|13p . in certain 
cases this form might be better suited for the particular problem. Suppose for example that the values of / 
are only available at lattice points, and as such performing the differentiation prescribed in Eq. (|13b[) cannot 
be carried out explicitly. For that, first a smooth multivariable interpolation would be necessary, but this 
intermediate (auxiliary) step is totally superfluous, as Eq. (fTl)]) yields the same result, thereby revealing the 
true usefulness of this formula. As to the solution of the Schrodinger equation, which is our primary object, 
one must remember that the values of / will indeed be available at the lattice points only, as according to 
Eq. |4b| 

fpikl = (E — U p ikl)ljjpikl, (17) 

and this results in the usual self-consistent equation, the quantum mechanical eigenvalue problem. 

Equation (|16j) is the main result of this subsection. In order to find the eigenstates and the energy eigen- 
values we need to say something about the remaining parameters, because they might affect the outcome. 
We will see indeed that they do. Therefore now we proceed with the analysis of this question. 



2.3. Free parameters and matrix representation 



In this subsection our main concern is the issue of free parameters 7, S and 7', that are left undetermined 
in the Schrodinger equation, Eq. (|16[) . This equation, in conjunction with Eq. (jTTJ) , is in fact a huge coupled 
system of linear difference equations, where the indices of lattice points run in the range 1, . . . ,n. These 
points make up a dense cubic grid inside the 4D hypercube. Lattice sites, that are on its surface, have 
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coordinates where at least one of p, i, k or I is equal to or n + 1. According to the boundary condition of 
Eq. (|4d[) , here ippiki = 0. Further, the physical wave function must also vanish outside the boundary because 
of the infinite potential barrier we have imposed. This we know from the analytical solution. We have already 
pointed out these before, but at this point it is worth repeating and having them in mind because of the 
following intriguing property of Eq. (|16j) . The fourth- neighbor operator, when acting on a function, takes 
into account neighboring sites that are two steps away along the coordinate axes from a given point. If we 
evaluate O^ippiki right next to the boundary, that is say p = 1, the result involves values of ip taken outside 
the boundary, in this case i/'-Mfcz- Function values taken at grid points with indices equal to — 1 or n + 2 
are fictitious. If we were to write down the equations that apply to them, new values were introduced that 
are located even farther from the boundary, leading to an infinite hierarchy of equations. We would like to 
emphasize that though it is very tempting to set these values to zero (because of the potential barrier), and 
thereby cutting the hierarchy, from numerical point of view this procedure is in principle not adequate. The 
true physical wave function does indeed vanish identically outside the domain, the non-zero fictitious values 
of ip, however, serve to determine the field inside consistently. The only condition we can make use of when 
solving the equations for the field inside is the boundary condition. 

In order to cut the infinite hierarchy and resolve this problem one can invoke Lagrange interpolation and 
express fictitious values with those that lie inside the domain. To retain full consistency in the sense that 
errors introduced in different ways are at most 0(h e ) a six-point interpolation is necessary, for example [19) 



.7 = 1 V/ 



ipj-i,ikl- (18) 



The outlined interpolation technique should only be applied when Eq. (116|) is considered right next to the 
boundary and even in this case the extra precision it provides compared to that when all fictitious values 
are set to zero by hand is presumably negligible, given the fact that ip vanishes on the boundary anyway. 
Therefore, hereafter we neglect this and take ippiki = everywhere outside the surface of the hypercube. It 
will turn out soon that luckily this step does not affect the required 0(h e ) accuracy of our calculations, at 
least as far as the eigenvalues are concerned. We note that these kind of difficulties arising at the boundary 
are well known in lower dimensions as well, not just in Numerov's method [9114) but also in the usual classical 
fourth-order and sixth-order methods [5]. 

According to what has been said so far, it is clear that among the difference operators CI4 is somewhat 
special. Moreover, as Eq. (|15cj) suggests, the prefactors S and S' cannot be made vanish at the same time, 
otherwise 7 = would hold too and we would be left with only the trivial solution of Eq. (fTT]) . So, at this 
point we have a degree of freedom how we distribute O4 on the two sides of Eq. (fTO)) . We find it convenient 
to set 5 = 0, thereby eliminating it from the left side. Then, knowing 7^0 must be, we can take without 
loss of generality 7 = 1 to obtain 

- ^ (12Di + n 2 + D 3 ) ippiki 

= U+U27'-ijni + (i-47'Jn 2 +7'n3-^on4)/ P iw, P ,i,k,i = i,...,n. (19) 

Now it is easy to see that all four difference operators appear indeed in the equation, neither of them 
can be completely eliminated, provided of course that we insist on the 0(h 6 ) accuracy. From this and from 
Eqs. © it follows that the method we have just developed is actually a fully discretized 89-point central 
difference scheme on the 4D grid. As mentioned in the introduction in Section [TJ in the 2D Schrodinger 
problem the full and partial discretization methods have already been applied for the study of anisotropic 
Coulombic potentials 7 9 10 . However, in dimensions greater than two the literature is not very helpful on 
this, to our knowledge no detailed calculations have been performed so far [20] . Our approach might be one 
of the first to fill up this gap. 

Setting values of two parameters turned out to be quite easy, but it is not the case with the last one and 
the rest of this subsection is devoted to this issue. In order to show how nontrivial the result is compared 
to 7 and <5, we shall anticipate it here 
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^ = 3840 * °-°° 59 - (20) 

Before we go into details of its proof, we believe it is worthwhile and instructive to examine how Eq. (|19|) 
would be modified if we followed the standard approach and derived it with an error of 0(h 2 ). This would 
then be the analogue of the most simple five-point scheme in two-dimensions [7 9 19J. It is not hard to see 
that in this case only the first term remains in the large parentheses. This corresponds to the usual case of 
an eigenvalue problem where there is only the identity operator on the right hand side. As opposed to this, 
had we wanted to reach an error of at most 0(h 4 ), it would have been sufficient to stop at the second term 
in Eq. (fT3b|) and to express that with difference operators. This result leads us to the following important 
conclusion: whenever the required precision reaches 0(h 4 ), or higher as in Eq. ([TO]) , the usual eigenvalue 
problem turns into a generalized one because even the right hand side, where the eigenvalue E shows up, 
contains difference operator (s). This observation characteristic to Numerov's method is well known in lower 
dimensions as well }6|9|10j . 

The difference equation (|19|) means in fact n coupled linear equations. In order to handle them together 
we form large column vectors, so-called stacks of dimension n 4 as 



= Vpiki, (21) 

where the one-to-one correspondence between lattice point indices and the vector index (/i = 1, . . . , n 4 ) is 
given by 

fi = p + n(i-l) + n 2 (k - 1) + n 3 (l - 1). (22) 



By this procedure any function v defined on the 4D grid can be mapped into a vector v € R" and vice versa. 
Equation (|22|) might remind one of the relation between indices of a Kronecker-product (direct-product) of 
matrices and those of the constituents. This is not surprising as v can be written as 

n 

v = v pik il(g>k®i®p, (23) 

p,i,k,l— 1 

where p, i, k and 1 are column vectors of size n whose elements differ from zero only at the pth, ith, 
kih and Zth position, respectively, where they all equal one0 Now that we know how to map multivariable 
functions to column vectors, it is obvious that linear operators acting on grid functions are isomorph to linear 
transformations of K" . As a result, the difference operators of Eq. © can also be mapped isomorphically 
to Mi € M n 4[R] matrices, where i = 1, ... ,4 and M„4[R] denotes the set of all n 4 -by-n 4 real matrices. 
Straightforward but lengthy calculations yield 

□ i < — ► Mi = A (g) E„ 3 + E n ® A ® E„2 + E„2 ® A ® E„ + E„ 3 ® A - 8E„4 (24) 

for the first-neighbor, 

□ 2 < — > M 2 = A ® A ® E„2 + A®E n g)A®E„ + A® E„ 2 ® A + E„ ® A ® A <g> E n 

+ F, n (g)A(g)~E n (g>A + E„2 <g> A ® A - 24E„4 (25) 

for the second-neighbor and 

□ 3 < — ► M 3 = A®A(g>A®E n + A®A®E n (g)A + A(g>E n (g>A(g)A + E ri (g)A(g)A(g>A - 32E„4 

(26) 

for the third-neighbor difference operator, respectively. Here E^ is the unit matrix of size k and A is the 
n-by-n tridiagonal matrix 



1 Note that two definitions are used commonly in the literature for the Kronecker-product of matrices, C = A ® B. We use 
the one that places the second matrix in the first, that is C is a large block matrix (sometimes called hypermatrix), where the 
ith block in the jth column is OjjB. 
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1 

1 1 



1 1 



(27) 



V 10 / 

Before we proceed with M4, we would like to call the attention to an interesting observation. Namely, the 
matrices above, if written as direct-products, exhibit similar "patterns" as the corresponding expressions 
in Eqs. ([6]). It is not so hard to see that the very simple tridiagonal form of A corresponds to that when 
values of ippiki are considered and added at first-neighbor sites in one given dimension. Then, by means of 
direct-product, linear combinations of any neighbors of a given site in any lower dimension subspace can be 
constructed. Coming back to the fourth-neighbor difference operator, this reasoning leads us to 



□ 4 < — ► M 4 = A<g)A<g>A<g)A - 16E„ 
+ A' ® E„3 + E n ® A' <g 



E, 



E r 



A' ® B r . 



E,. 



A' - 8E„4 



(28) 



where A' is n-by-n and corresponds to combination of second-neighbors along a given coordinate axis. The 
second line on the right hand side is the matrix representation of the last term in Eq. (|6d|). Comparison 
with Eqs. I|24p and (|2T|) suggests A' should be pentadiagonal 



1 

1 

1 1 



10 



V 



1 



(29) 



Here, in the first and last two rows there is only one nonzero element. This is because along any given 
coordinate axis the first two and last two grid points each have only one such second-neighbor that is also 
inside the domain. Right on the boundary ip P iki = because of the boundary condition, whereas, as pointed 
out in the beginning of this subsection, fictitious values of ippiki were taken to be zero by hand. It was also 
emphasized there that this procedure does not affect required precision, which we will prove exactly at the 
end of this subsection. Now we make another presumably very accurate approximation that does not affect 
precision either: instead of Eq. (j29|) we take 



A' 



2E K 



V 



1 1 

1 

1 1 

10 
1 -1 



\ 



(30) 



that is A' u = A' nn = —1, otherwise A' is the same as Eq. (|2"9"|) . The sparsity of Mj can be seen in Fig. [TJ 

Having obtained these important matrices we are now in a position to vectorize Eq. (|19p . Taking into 
account Eq. (fT7|) as well we finally arrive to the 0(h 6 ) matrix form of the discretized four-dimensional 
Schodinger equation 

(/i- 2 M + Ndiag(U))i/> = ENtl>. (31) 
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Fig. 1. (Color online) Sparsity patterns of Mi, M2, M3 and M4 (from left to right) at a very low resolution of n = 3. Their 
total size is n 4 = 81. 



Here, diag(U) is a diagonal matrix composed of the vector U and the matrices M and N, which include all 
information supplied by the boundary condition, arc 

1 



M 



30 



(12MX + M2 + M3), 



and 



N = E.„4 



12 7 ' - — ) Mi + ( — - 4V) M 2 + 7'M 3 - — M 4 . 
' 30 / V 36 ' ' 240 



(32) 
(33) 



Apropos of Eq. (fl9|) we have already mentioned that from linear algebraic point of view Eq. (|3lj) constitutes 
a generalized eigenvalue problem because N ^ E„4. In what follows, we collect its basic mathematical 
properties as well as the necessary requirements it must satisfy in order to best represent the continuous 
Schrodinger equation, Eq. ((4]) . We shall see that these will lead us to the very specific choice of 7' anticipated 
in Eq. ([2U]). 

(i) Properties of the Kronecker-product and the symmetry of A imply that M, (i = 1, . . . , 4) as well as 
M and N are all symmetric matrices. 

(ii) As shown in Appendix [Al the matrices are all negative definite. From this, using the properties of 
definitcE] matrices, it follows that M is positive definite. 

(iii) Let 7' be such that N is non-singular. Then multiplying Eq. (f3"Tj) by N _1 one ends up with the usual 
form of a quantum mechanical eigenvalue problem where the kinetic energy operator — A is represented 
by /i _2 N _1 M. It is also shown in Appendix [Al that M and N commute: [M, N] = 0. This is a very 
important and central result for it assures that the matrix Hamiltonian (/i~ 2 N _1 M + diag(U)) is 
symmetric leading to real energy eigenvalues. 

According to (iii), /i~ 2 N -1 M represents the negative Laplace operator. It is known from functional 
analysis that —A is a nonnegative and non-bounded linear operator, thus 7' must be chosen so that 
N _1 M is positive definite [11] . Again, using the properties of definite matrices and taking advantage 
of the fact that the terms of this product commute, beyond the requirement for N to be non-singular 
it must also be positive definite. In addition to that, asymptotically, as n — > 00, the spectrum of 
/i _2 N _1 M must approach that of —A supplied by the boundary condition 



(iv) 



lim ^(n+lfN-'M) 



1,2. 



(34) 



Note that the right hand side is nothing else but the set of allowed discrete energy levels of the 

four-dimensional "particle in a cubic box" problem. 
In Appendix [Al we present detailed derivations of the statements in (ii) and (iii). Also, we give analytical 
proofs of the followings. The necessary requirement in (iv) for N to be positive definite can be satisfied if 7' 
is taken from the interval 

23 



7 



< 



3840' 



(35) 



We say that a matrix is definite if it is either positive or negative definite. 
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Besides, the same direct polynomial structure of M and N imply that all n 4 eigenvalues of N _1 M are of 
the form 0(k)/A(k), where 0(k) and A(k) are eigenvalues of M and N, respectively. Explicit expressions for 
these quantities are given in Eqs. (|A.13|) and (|A.26[) . Here, k = (k±, &2, ka, fej,) is a composite vector index 
with components taking integer values in the range 1, . . . , n. With these results at hand we can perform 
asymptotic expansion for large n, and by Eq. (|A.37j) we find 



, n .U, <k)-6(k) 7 ' 



(n+1) AM = e(k) ' 1 ~ -» +•■■ )' (:!(,! 



where 



4 

-2 x 

»=i 



e(k)=7r 2 Vfc| ) ^ = 1,2,... (37) 



is the energy level in Eq. (fM)) . This result forms the cornerstone of our 0(h 6 ) theory, for it assures that 
(i) as n goes to infinity, the spectrum of the discrete Laplace operator evolves into that of —A, (ii) the 
leading error is at most 0(n -6 ). We see now finally that the approximations we applied by setting fictitious 
function values to zero by hand and using Eq. (|30p instead of Eq. ((29)) are indeed that accurate as were 
expected. They did not affect the required precision. This is a very satisfactory situation because it suggests 
that higher order calculations could be performed in exactly the same manner. 

Tuning 7' within the range of Eq. (|35[) the order of error cannot be reduced any further. However, we can 
still try to minimize its prefactor. According to Eqs. (|A.38|) and (|A.39|) this can be achieved successfully by 
taking the largest value allowed 

->' - SS- (38) 

and this is exactly what we have anticipated in Eq. (|20p . It might be of interest that with this choice N 
becomes asymptotically singular, suggesting numerical implementation should avoid manipulation with the 
inverse. For more details please, see Subsection IA.4I of Appendix lAl 

Now that we have successfully derived the matrix Schrodinger equation and built up the relevant matrices, 
in the last subsection we give a very brief discussion on the memory demand of the storage of these matrices. 



2.4. Sparsity properties, memory consumption 

Figure[T]shows the sparsity patterns of the matrices M, , that are the basic building rocks of the Schrodinger 
equation, Eq. (|31j) . As for the resolution, for illustration purposes we used a very low value of n — 3, which 
means that the total size of the matrices is n 4 = 81. Remarkably, we found that due to the direct-product 
structure of these matrices the number of nonzero elements and therefore the density or the sparsity can 
be calculated exactly for arbitrary n. Since the number of nonzero elements in A equals N(A) = 2(n — 1) 
it follows form Eq. (|24| that for Mi it is 8(n — l)n 3 + n 4 . Similarly, for M 2 , M 3 and M 4 the number of 
nonzero elements are exactly 24(n — l) 2 + n 4 , 32(n — l) 3 n + n 4 and 16(n — l) 4 + 8(n— 2)n 3 + n 4 , respectively. 
Knowing these values it is now straightforward to obtain respective results for M and N and they read 

N(M) = 8(n - l)n 3 + 24(n - 1)V + 32(n - l) 3 n + n 4 w 65n 4 , (39a) 
N(N) = 8(n - l)n 3 + 24(n - 1) V + 32(n - ifn + 16(n - l) 4 + 8(n - 2)n 3 + n 4 w 89n 4 , (39b) 

where the approximate values are the leading terms for large n. It also turns out that the matrix in the 
parentheses on the left hand side of Eq. (f3Tj) has exactly the same structure as N (from the sparsity point 
of view), which means that for large n one has to allocate memory for roughly 178n 4 nonzero real matrix 
elements in order to simply store the equation as a whole. For example at a modest resolution of n = 20 and 
with double precision arithmetic on a 64 bit computer the memory consumption is roughly 0.37 gigabytes. 

In this Section we have established the algebra of the two-particle problem, it is now time to proceed 
with the group theoretical aspects. These, as expected, result in further substantial simplifications that will 
facilitate numerical implementation. 
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3. Group theoretical aspects: classification of energy levels 

In the previous section we have developed a new finite difference method for the numerical solution of the 
two-particle stationary Schrodinger equation. The most noticeable feature of the theory is that it provides 
the energy levels very accurately. The physical problem we are considering in this paper is the interaction 
of two identical particles in two-dimensions and with closed boundary condition. A closed boundary in our 
case means that both particles are confined in an infinite square well, which leads to the Dirichlet condition 
in Eq. (|4d|) . For these circumstances we have derived a matrix Schrodinger equation in Eq. (|3"T|) . We showed 
analytically that as the step size h goes to zero it converges to the continuous Schrodinger equation with 
a local truncation error of 0(h 6 ). The problem can be thought of as an eigenvalue problem of an abstract 
four-dimensional particle in a 4D cubic box. From this and from the details of derivation one can see that 
the algorithm can be applied in any dimension, or, if necessary, the precision can be refined further as well. 

Despite of the favorable high precision there are, however, difficulties too. When it comes to numerical 
solution, we face a generalized eigenvalue problem of dimension n 4 . As we saw in Subsection 12.41 due to the 
multiple direct-product structure even a modest resolution in n results in huge matrices that, though being 
sparse, still require lots of memory. Also, what is more important, the classification of energy levels and 
wave functions is unknown at this stage. It might happen that we unnecessarily calculate degenerate levels 
many times and cannot reach higher energy excited states that might be of interest. The good thing is that 
these problems can be entirely overcome as there exists a powerful theoretical tool, group theory, that was 
developed exactly for these needs. 

Before jumping into details we would like to emphasize already at this point that the group theoretical 
analysis we are about to present here is not by all means necessary in order to solve Eq. ([3T|) . With a suffi- 
ciently powerful workstation the diagonalization at this stage is in principle straightforward 20J. However, 
exploiting all the hidden information that lies in the symmetry of Eq. (|31j) is very instructive an far reaching 
as we shall see right away. 

3.1. Group of the Schrodinger equation: general remarks 

The very first step in the development of representation theory is the exploration of the maximal invariancc 
group Q of the Hamiltonian H. In quantum mechanical description of a system Q consists of all unitary 
operators of the Hilbert space that commute with H. In many cases it turns out that these groups are, 
in group theoretical sense, linear Lie groups and it is the generators of the associate real Lie algebra that 
represent the physical quantities that are constants of motion. As the generators are Hermitian operators that 
also commute with H , they are indeed the usual conserved quantities of the system. It is a very interesting 
observation that the exploration of the maximal invariance group can be quite complicated sometimes 
[2"T] . Perhaps the two most famous problems of quantum theory are the hydrogen atom (also known as the 
quantum mechanical Kepler problem) and the isotropic harmonic oscillator [22123] . As a matter of fact, both 
problems can be considered in any dimension n. These problems have an obvious geometrical symmetry: 
due to rotational symmetry of the central force, Q, whatever it may be, must contain the n-dimensional 
orthogonal group 0(n, R) as a subgroup!^] This immediately explains degeneracy of the energy levels of 
both systems in the angular hypermomenta rrii, i — l,...,n — 2, but cannot account for the degeneracy 
in the total angular momentum I. Detailed investigations show that behind the substantial "accidental" 
degeneracies there are larger symmetry groups: in case of hydrogen atom it is 0(n + 1), whereas for the 
isotropic oscillator SU(n) is found 22 24 25 26J. Even the respective classical Hamiltonians possess these 
enlarged invariance groups. The extra symmetries spanning the larger groups, that were initially overlooked, 
are the so-called hidden or dynamical symmetries that cannot be associated with apparent symmetries of the 
configuration space |21l27j . They are related to the fact that the Schrodinger equation can be separated in 
other coordinate systems as well as in spherical polars. The additive generators of the respective Lie algebras 



3 Hereafter we drop extra notation R, because we only consider real matrices throughout the paper and there will be no 
confusion by simply writing O(n). 
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are in turn the Runge-Lenz vector operator and a symmetric tensor operator in case of the hydrogen atom 
and the oscillator, respectively [22l28j . 

Our main concern in this section is the exploration of the symmetry group and development of represen- 
tation theory of the interacting two-particle problem under focus. This, as we have already noted before, 
can be looked at as a one-particle problem in a four-dimensional space where the particle is trapped in an 
impenetrable cubic box. In addition to that there is a well defined potential inside the domain too. The 
problem has an obvious geometrical symmetry due to the cubic structure of the cavity. Therefore Q will 
certainly contain discrete symmetry operations related to the symmetry of the configuration space and the 
external potential. The more interesting question, however, addresses the possible existence of dynamical 
symmetries. In order to give a full classification scheme of the levels and wave functions it is necessary to 
incorporate these symmetries as well, if there are any. It turns out that the Coulomb interaction we are 
about to apply excludes the possibility of such hidden symmetries and we are left with a discrete symmetry 
group consisting of geometrical symmetries only. In this specific case the unitary operators of the Hilbert 
space that build up Q are the so-called scalar transformation operators [29]. We give detailed analysis of 
them in Subsection l3.3l These operators are known to be isomorphic to the geometrical symmetry operations 
of the system, thus Q can equally be considered as the group of abstract coordinate transformations of the 
four-dimensional space. In what follows we find it more convenient to use this latter designation. 

Without interaction the problem, though elementary from quantum mechanical point of view, from the 
viewpoint of group theory becomes much more elaborate. The lack of interaction leads to substantially larger 
number of discrete symmetries and more interestingly to the appearance of dynamical symmetries as well. 
This issue will be partly explored in Section |U 

3.2. Group of the Schrodinger equation of the confined interacting particles 

According to what has been said in the second paragraph of the previous subsection, this group can be 
considered as consisting of those four-dimensional coordinate transformations T that leave H invariant. To 
put this requirement in a mathematical formula recall that H is now 

fl w = -E£+^ ( 4 °) 

l—\ L 

where the potential is given by Eq. ((3]). As is known from quantum mechanics, the n-dimensional Laplace 
operator is invariant under pure rotations of the orthogonal group O(n). In our case n = 4 of course. Let T 
be a coordinate transformation that introduces a new set of mutually orthogonal axes (x^ , x' 2 , x' 3 , x\ ), but 
which leaves the origin unmoved. Now there exists a unique matrix R £ 0(4), with which the coordinates 
of a fixed point with respect to this new system can be expressed as x' = Rx, with x being the coordinates 
with respect to the original frame (x±, X2,xs, x<l). The transformation T is then said to leave the Hamiltonian 
invariant if 

iJ(Rx) = H(x). (41) 

The group of abstract transformations T is isomorphic to the group of matrices R, hence we will refer to 
coordinate transformations as matrices and Q will hereafter be considered as the group of the corresponding 
matrices. 

Due to the presence of the potential and the strict boundary condition our system is not homogeneous, 
pure translations do not leave H invariant. Further, if R S Q is a symmetry transformation, then (i) it must 
be a symmetry of the four-dimensional hypercube and (ii) satisfy 

£7(Rx) = U(x). (42) 

As to the symmetries of the domain, the n-dimcnsional hypercube has altogether 2 n n\ symmetry transfor- 
mations, all of which can be faithfully represented by n-by-n orthogonal matrices having only one nonzero 
element, +1 or —1 in each row and column. Such matrices are called signed permutation matrices. It follows 
there are 384 such matrices in 4D, these constitute the so-called four-dimensional cubic group O4 [3D]. To 
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single out those that make up Q we have to specify the potential too. The external field the particles are 
subjected to is chosen as 



U(r) = ^ 2 (x 2 +y 2 ), 



(43) 



where uj > is a dimensionless parameter responsible for its scale. We use quadratic potential because for 
this the noninteracting problem has exact analytical solution 



E = cj (v kl + v k2 + Uk 3 + Vk 4 + 2) , 



(44a) 
(44b) 



Here k = (fci, fe, fe, k±) and each quantum number takes integer values: foj = 0, 1, 2, . . . . The wave functions 
<t>k and the first sixteen z/^'s are shown in Eqs. (jB.llj) . (1B.15|) and in Table |B~T1 of Appendix IB] respectively. 
Besides it is not singular as opposed to atomic potentials (the 2D hydrogen atom of Ref. [3T] for example) , 
which is also convenient from numerical point of view [2 18 32 . With this choice we find that the first two 
terms of U together are invariant under all rotations of 0(4). In contrast to this, the pair interaction V 
constitutes a much stronger constraint: it is invariant under such orthogonal transformations for which 



ri - r 2 



where 



(l 



Q = 



x T Qx 



0-10 



10-1 
-10 10 







-1 1 



(45) 



(46) 



Now, a transformation R can only be a member of Q if [R, Q] = 0, which in turn restricts the 384 potential 
candidates to a subgroup of order 32. This can be verified either by direct calculation or by the following 
alternative method. Consider the group A of order 4 with members 



Ri 



l 



\ 



V 



1 



/ 



Ro 



1 



V 



1 



Rs 



/ 



v 



i 



i 



A 



Rt 



V 



1 



(47) 



/ 



Only the nonzero elements are shown for better readability. These matrices correspond to permutations of 
coordinates of the particles. For example R 2 permutes the two particles (xi,^) = (^i 7 yi) *■ — ► (%2,1/2) = 
(x3,X4), whereas R3 and R4 permutes only the x and y coordinates, respectively. On the other hand, the 
symmetries of the confining two-dimensional square well form another group B of order 8 with elements 



(l 



Ri 



Rr; 



\ 



R 



2 — 



1 



Ri, 



Ri; 



R2 



Ri 



R 7 



-1 



R 



4 — 



1 



-1 



Rr 



-R4 



(48) 



We can observe the block diagonal structure of these matrices with the same 2-by-2 blocks repeated twice 
in the diagonal. This reflects the fact that the same two-dimensional transformation must be performed on 
both particles. We note that B is essentially a four-dimensional faithful representation of the group of point 
symmetries of a square, denoted commonly by C<± v in solid state literature. Every member of A and B obeys 
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Eq. (141 [1 and so do their products. This results in the rather plausible observation that the group of the 
Schrodinger equation is expressible as 

Q = {R Q R 6 I R„ e A, K b € B} . (49) 

It is indeed a group of order 32, where A and B are trivial subgroups. Closer inspection shows that A is 
an invariant Abelian subgroup of Q and isomorphic to C|, with C2 being the cyclic group of order 2. This 
feature, in conjunction with the defining equation (|49[) and the fact that A and B have only the identity 
element in common lead to an important result: the full invariance group of the Hamiltonian has a semi-direct 
product structure [29j 

Q = A©B ~ Ci®C 4v . (50) 

This will come in handy when we proceed with the determination of the irreducible representations of g. 
Technical details concerning this can be found in Appendix [Cl 



3.3. Scalar transformation operators 



Now that we have successfully explored the group of the Schrodinger equation we go on and construct the 
scalar transformation operators. From group representation theory it is known that these linear operators 
are absolutely necessary for the explicit determination of symmetry adapted basis functions [29] . Let R € g 
a symmetry transformation of the Hamiltonian, then the scalar transformation operator -P(R) is defined by 



P(R)w(x) = w(R -1 x), 



(51) 



where w(x) is any multivariable function. Recall now that in the numerical analysis the values of a function 
are only available at 4D grid points, see Eq. ([5|). Thus a function is essentially nothing else than a huge 
set of real numbers conveniently mapped to column vectors by means of Eq. (|23|) . Accordingly, every linear 
operator acting on this space stands in one-to-one correspondence with a matrix of M n 4 [M] , for which we 
have already seen examples in Subsection 12.31 by the construction of the matrices M^. As to the scalar 
transformation operators we find 

{P(R)} p „ = 1, (52) 

where 



H = p' + n(i' - 1) + n 2 (k' - 1) + n 3 (l' 
u=p + n(i-l) + n 2 (k - 1) + n 3 (l - 1), 

and the relation between indices is given by 



1) 



(53) 
(54) 



n + 1 



Ei. - R) 



1 




P 


1 




i 




+ R 




1 




k 


1 




I 



(55) 



Running the indices p,i,k and I in the range 1, ...,n the whole matrix can be constructed. Further, by 
Eq. (|52p these have exactly one nonzero element in each row and column, they are necessarily orthogonal 
and describe permutations. This must be so as P(R)v is by definition the same scalar field as v, but now it 
is looked at from within the transformed coordinate system introduced by R. In other words we can say that 
every such matrix is a unitary linear transformation of K™ . Similarly, the P(R)'s are unitary operators of 
the Hilbert space [29] . In complete analogy with Eq. |23|) they can be expressed with multiple direct-products 
too 



P(R)= (l'ol T )®(k'ok T ) 

p,i,k,l— 1 

where o stands for the dyadic product (outer product). 



(p' 



o p 



(56) 
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One of the most important theorems of group representation theory, at least from quantum mechanical 
point of view, is that the scalar transformation operators form a group isomorphic to the group Q of 
transformations R [29]. Further, as Eqs. (|4Tj) and (|5Tj) suggest, H and P(R) commute for all R £ Q. As we 
pointed out in Subsection l3.ll in quantum mechanics the group of such unitary operators is usually identified 
as the invariance group of the system. Now it is not surprising that the matrices representing these operators 
in the finite difference method obey a very same equation 

/ l - 2 N- 1 M + diag(U),P(R)j = 0. (57) 

Beyond this we can also prove that P(R) not only commutes with the full Hamiltonian, but does it also 
with each term separately. The fact it is true for the potential is almost trivial from Eq. (|42|) . while the proof 
for M and N is as follows. First remember that the basic building blocks of the finite difference method are 
the difference operators given by Eqs. ©. Every such operator can be expressed concisely as 

□i«(x) = X>(x + d) (58) 

d 

where the sum is over all jth neighbors d that are at a well defined distance from the given lattice point x. Let 
their total number be denoted by gj . In the four-dimensional grid we use for numerics this is equivalent to the 
constraint d 2 = h 2 (a 2 + /3 2 +7 2 + <5 2 ) = const., with a, . . . , S being integer coordinates of d. This expression 
is clearly invariant under any permutation of coordinates with sign changes included. In other words it is 
invariant under all 384 orthogonal transformations of O4 and as such for all R G Q as well. This leads 
immediately to the fact that the neighborhood is completely symmetric: {di, . . . , d gj } = {Rdi, . . . , Rd gj }. 
This is to be interpreted as an equation of sets where the order of elements is irrelevant. With these findings 
and the definition in Eq. (|51ll it is now easy to verify that 



[□i,P(R)]=0. (59) 

Taking into account Eqs. (152")) and (|3"3")l and the fact that the commutator is a real bilinear function, com- 
mutativity of the respective matrices is found 

P(R)MP(R)" 1 = M, (60) 
P(R)NP(R)^ 1 = N. (61) 

This fundamental result, which we shall make use of shortly in Wigner-Eckart theorem [29j in the following 
subsection, shows that beyond the potential and the full Hamiltonian these are also irreducible tensor 
operators (matrices) of the completely symmetric irreducible representation T 11 of Q. Representations and 
their labeling convention can be found in Appendix [Cl In particular, T 11 is defined by Eq. (|C.15jl . 

The final step in the development of group theoretical background is the determination of projection 
operators and the construction of a set of symmetry adapted basis functions. 



3.4. 



Projection operators, basis functions of irreducible representations and the low energy subspace 



In the previous subsection we have obtained matrix representations of the scalar transformation operators. 
Also, in Appendix [C] we have elaborated all irreducible representations T qp of the group of the Schrodinger 
equation. With this knowledge we can now construct a special orthogonal basis of the Hilbert space, every 
member of which is a basis function (in group theoretical sense) transforming as some row of some irreducible 
representation of Q . This is achieved by means of the projection operator method. Let's see how! 

The projection operators are defined by [29) 

= (<Vs) E r 9P (R) y P(R), (62) 

Res 

where g — 32 is the order of the group and d qp is the dimension of T qp . In general, these objects as well as the 
scalar transformation operators are linear operators acting on functions of the infinite dimensional Hilbert 
space H. In the finite difference method, however, since Ti. = M" they become large but finite dimensional 
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sparse matrices. According to developments of representation theory if v is a vector of this space such that 
Wj P = P'^v, j = 1, . . . , dq P are all nonzero, then they form a basis for T qp 

P(R)wf =gr"(R) y wf. (63) 
i=l 

Furthermore, the powerful Wigner-Eckart theorem says that if S is any irreducible tensor operator trans- 
forming as the completely symmetric irreducible representation T 11 , that is commutes with all P(R), then 
a remarkable simplification of the matrix elements occur 

(wf , Swf ) = 6 jf S qp , qlpl (wf , Swf ) , (64) 

and the scalar product on right hand side does not depend on the row index j. As to the real scalar product 
it is defined as usual: (a, b) — J2i a i^i- ^ n the previous subsection we showed that each term appearing in 
the matrix Schrodinger equation exhibits the same transformation property as S above. Now it follows that 
the eigenvalue equation of dimension n 4 can be transformed into an equivalent block diagonal form where 
each block belongs to some row of some irreducible representation of Q. Also, it is shown in Appendix [Cl that 
Q possesses 8 one-dimensional and 6 two-dimensional representations. In the block diagonal decomposition 
there will be therefore 8 + 2 x 6 = 20 blocks among which only 14 are different. The solution of the eigenvalue 
problem is then equivalent to the solution of each block separately. 

Let us now define a basis (in linear algebraic sense) of the Hilbert space. This will serve as input from 
which the symmetry adapted basis will be projected out. It reads 

(in \ ~ 5 

nE^i^''!') 4>k 1 (xi tP )(f>k 2 (x2,i)(l)k 3 (x3,k)(l>k i ( x ^l)^ (65) 
i=i P =i / 

where the quantum numbers ki take integer values in the range 0, ... ,m — 1 < n. Again, as in Eq. (|44b[) , 
the functions <f>k{x) are the exact eigenfunctions of the confined one-dimensional harmonic oscillator. Their 
analytical form is given in Appendix [Bl in Eqs. (|B.11[) and (|B.15[) . From these grid values v(k) is obtained 
by Eq. (f2"3")l . If m — 1 = n they constitute an orthonormal basis of the Hilbert space. On the other hand, if 
m — 1 < n, they generate a subspace Tt m C TL of dimension to 4 , that is physically of most interest. In what 
follows we solve the discrete Schrodinger equation in this low energy subspace. This we can do because we 
are primarily interested in the ground state and some higher energy excited states only and for this purpose 
it is completely sufficient to restrict our analysis to the low energy subspace. Of course, with sufficiently 
powerful computing facilities one can increase to and thereby expand the computational subspace. 
In the basis of v(k) the equation looks like 



E 

k' 



r(k), (V 2 M + Ndiag(U)) v(k')) V(k') = E ^ (v(k), Nv(k')) V>( k ')- (66) 



This is actually the kth row of the Schrodinger equation for the unknown coefficients ip(k) and energy E. 
The corresponding state is then obtained as 

</> = 5>(k)v(k). (67) 

k 

Due to the scalar product the matrices of Eq. (|66p are not sparse anymore and still quite large. Also, the level 
classification is still unsolved at this stage, hence this is the appropriate point to appeal to the symmetry 
adapted basis {wf (k)} instead of {v(k)}. This is obtained by applying the orthogonal projections 

wf(k)=Pfv(k). (68) 

From what has been said so far in this subsection it is clear that this newly generated vector is necessarily 
contained in TL and, if not zero, transforms as the jth row of T qp . However, it is not obvious that (i) it is 
also a member of TL m and (ii) the maximal linearly independent set of these is orthogonal and span 7i m . 
These important properties are proved in detail in Appendix [Dl 
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Consider now a representation T qp . In Appendix iDl we give an algorithm to generate a maximal set of 
independent vectors that all belong to the jth row of this particular representation. These will be denoted 
by w^ p (k s ), s = 1, . . . , r qp . With all these the Schrodinger equation in this symmetry channel is found to be 

r qp r qp 

J2 (wf (k s ), (V 2 M + Ndiag(U)) wf (k t )) A = (wf (k s ), Nwf (k t )) 8=1,..., r™, 

t=i t=i 

(69) 

and the actual wave function of this symmetry is 

Vf = $> t wf (k t ). (70) 

t=i 

Equation (|6"9"|) is one block out of the twenty of the block diagonal form of Eq. (|6"6")) . As the scalar products 
do not depend on j, energy levels obtained from this will be exactly d gp -fold degenerate. In the numerical 
implementation we used n = 30 and m = 8, see Sections Q] and O This means Ti and TC m are of dimensions 
n 4 = 810000 and m 4 = 4096, respectively. This latter dimension decomposes into smaller blocks of sizes r qp 
in our case, each labelled by a certain representation. These sizes can be be determined either theoretically 
by means of of Appendix [D] or numerically. Our computations showed that these are at most a few hundred 
(4096/20 « 200 in average), which are very easy to handle numerically. Precise values are shown in the 
second row of Table [TJ 

3.5. Permutation symmetry of the two-particle wave function: Fermi or Bose statistics 

So far we were not concerned about the nature of particles we are dealing with. That is we were not 
interested in whether they obey Fermi or Bose statistics. The only constraint we ordered during the problem 
formulation in Subsection 12.11 was that they are completely indistinguishable. From this it follows that the 
Schrodinger equation has both symmetric and antisymmetric solutions. Also, as is well known from quantum 
mechanics, fermions have a total wave function (including spin), that changes sign whenever two particles 
are interchanged. In contrast to this the total wave function of bosons is completely symmetric. However, 
since in this paper spin related effects are not considered, we omitted the factorizable spin state in the very 
beginning and concentrated on the real space component only. Having developed the group representations 
its permutation symmetry is now easy to check as follows. 

The group element in Q describing permutation of particles is R2 € A, see Eq. (14"T1) . The scalar transfor- 
mation matrix for this operation is in turn 

n 

P(R 2 )= (i°l T ) ® (p°k T ) ® (loi T ) ® (kop T ) . (71) 

p,i,k : l— 1 

Knowing the irreducible representations of Q (from Appendix [U]) , one can easily verify that 

I^(R 2 ) = {- E2 ' i^ = 2andp = l,...,4 ) 
[+E dw , otherwise, 

where E^, as before, denotes the unit matrix of dimension k. These results in conjunction with Eqs. (|63p 
and ([70| show that all solutions of the discretized Schrodinger equation satisfy the transformation property 

P(R 2 )Vf = T^f • (73) 

From this we see that it is the representations that unambiguously distinguish between solutions with 
different permutation symmetry: eigenvectors transforming as some row of T 2p are all antisymmetric and 
span the antisymmetric subspace of the total two-particle Hilbert space. On the other hand, all other 
solutions that belong to the remaining ten irreducible representations are symmetric and span the symmetric 
complementary subspace. Note that each of these subspaces can be physically important: for example if the 
two particles under consideration are fermions of spin one-half (electrons, protons, neutrons, . . . ), then the 
spin singlet and triplet states are accompanied by symmetric and antisymmetric wave functions, respectively. 
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4. Numerical results for noninteracting particles 



This and the next sections are devoted to the demonstration of numerical results that are based on the 
theory presented so far. Here the noninteracting problem is concerned. In this case the quantum mechanical 
problem of the two identical particles can be traced back to a single-particle problem and analytic solution is 
easy to find. Though the lack of interaction turns the problem into exactly solvable, in group theoretical sense 
it has some peculiarities. Namely, the symmetry group of the system is enlarged considerably and contains 
dynamical symmetries as well. These and other geometrical symmetries will be broken by the interaction. 



4.1. Enlarged symmetry group: dynamical symmetries 



Numerical solution of Eq. (|69|) in each symmetry channel provides the level structure and eigenstates of 
the system. For noninteracting particles these are already known analytically in Eqs. ()44|) . Exact energy 
eigenvalues E/lo are shown in the first column of Table [T] to six digits of precision for parameter values 
of w 2 /2 = 500 and 6=1. They were calculated from Eq. (|44a[) and Table |B~T1 of Appendix IBl Numbers 
in round brackets indicate the degrees of degeneracy of each level. At this point we would like to call the 
attention to the fact that because of the four-dimensional cubic structure of the boundary condition the 
eigenvalues themselves as well as their degeneracies differ from those of the usual 4D isotropic oscillator. 
There the energy would be E — uj(n + 2), with n a nonnegative integer, and the degeneracy of the nth level 
is [24T25] 

fn + 3\ 

g(n) = ■ ( 74 ) 

The highly degenerate nature of the levels is in complete accordance with SU(4) symmetry of the problem 
[22j . In our case, however, the confining potential well breaks down full 0(4) rotational symmetry and leads 
to the finite four-dimensional cubic group O4. This partly explains for example why the degeneracy of the 
third level is not 10, but splits into 6 + 4. Similar splittings can be observed further down the first column of 
Table [TJ In solid state terminology this is nothing else then a crystal field splitting, where the perturbation 
in our case playing the role of the crystal field is the infinite potential barrier. 

A more interesting observation is, however, that the representation theory of O4 cannot explain completely 
all degeneracies found in the level structure, because the dimensions of its irreducible representations are 1, 
2, 3, 4, 6 and 8 [30]. This means, again looking at Table HJ that the twelve times degenerate energy levels 
either exhibit accidental degeneracies or beyond geometrical symmetries there are dynamical symmetries 
as well and the invariance group is actually larger than O4. Given the simplicity of Eq. (|44ap and the fact 
that since the z/„'s are related to the zeros of the confluent hypergeometric function they are definitely not 
integers, one might have the impression that accidental degeneracies are very rare or more likely do not occur 
at all. In the latter case simple combinatorial reasoning shows that the degree of degeneracy of a particular 
level can only be 1, 4, 6, 12 or 24. Following Ref. [33J we found indeed that extra degeneracies are related 
to the existence of a larger symmetry group and it is the irreducible representations of this group that are 
associated with the energy eigenvalues. This group will be explored next. 

The Hamiltonian without interaction can be written as H = ^ i Hi, where 

d 2 1 

Hi = ~dx 2 + 4 UJX i' i=l,...,4, (75) 

and as is known 

[H i ,H j }=0. (76) 

This signals it is not just the total energy that is conserved but also each component separately. In other 
words, due to the lack of interaction the one-dimensional sub-oscillators do not exchange energy quanta 
during the motion and this is essentially the reason for that the Schrodingcr equation can be separated in 
Descartes coordinates. Remember that the Schrodinger equation is supplemented with a boundary condition 
of Eq. (I4djl that prevents separation in spherical or cylindrical polars, thus 4D angular momentum is not 
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Table 1 

Classification of the first few energy eigenvalues E/lu of the two noninteracting particles in the confined harmonic potential. 
Scaling factor u) and the half box size 6 were chosen as w 2 /2 = 500 and 6 = 1. First column shows the exact theoretical 
results to six digits of precision calculated from Eq. I|44a| l and Table IB. II of Appendix [B] Numbers in parentheses indicate 
the degrees of degeneracy of each level. In the second column numerical data are given that were obtained from Eq. (I69D for 
m = 8 (the dimension of the low energy subspace is therefore m 4 = 4096) and for a resolution of n = 30. The step size is 
h = 26/(n + 1) = 0.0645. Further columns show the level distribution across representations of Q, whereas the last one accounts 
for the total degeneracy of the given level. Note that representations written in boldface are two-dimensional. The numbers 
in the second row mean the total number of independent vectors in the given representation. These we denoted by r qp in the 
text. Thus J2 d qp r" p = 4096. 
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a constant of motion. Neither are the offdiagonal components of the symmetric tensor operator found to 
be responsible for the extra degeneracies of the unconfmed oscillator |22|23j . In this subsection we refer to 
the invariance group Q' of the noninteracting Hamiltonian as the group consisting of all unitary operators 
of the Hilbert space that commute with H. We know so far there is a discrete subgroup, call it O2, that 
is related to geometrical symmetries: it involves the scalar transformation operators P(R), where R G O4. 
Commutativity of oscillator quanta with the total Hamiltonian suggests that beyond this there is a four 
parameter continuous subgroup as well consisting of 

[/(A 1 ,A 2 ,A3 ! A 4 ) = e i ^ Aj ^. (77) 

Due to commutativity of the generators in Eq. (|76|) the associated real Lie algebra is not semi-simple and 
the representation theory of the problem becomes quite cumbersome. Nonetheless, there exists a discrete 
group sufficient for our purposes. Let us define certain linear operators of the Hilbert space [33] 

A -(/>k(x) = sign(+j/ fel + v k2 + u k3 + ffcJV'k(x), (78a) 

Ai-0k(x) = sign(+i/ fel + v k2 - Vk 3 - ffcJV>k(x), (78b) 

A 2 -0k(x) = sign(+^ fel - v k2 + v k3 - ^ fc4 )?/> k (x), (78c) 

A 3 "0k(x) = sign(-^ fel + v k2 + Uk 3 - ^ 4 )-0k(x), (78d) 

where V'k(x) is the exact eigenstate of Eq. (|44bj) . These functions form an orthonormal basis of the Hilbert 
space, therefore Eqs. I|78j) in conjunction with linearity define the action of Aj on each function of the space. 
It is easy to see that since Vi are positive An is the identity operator /. Its definition reminds of the action of 
the Hamiltonian, albeit here it is not the energy that multiplies the eigenstate but its sign only. In a similar 
manner the other operators can be put in correspondence with linear combinations of Hi. Equations (|78[) 
lead to the following properties: (i) [Aj, fZ] = 0, (ii) [Aj,Aj] = 0, (iii) Af = I and (iiii) Aj are unitary 
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and self-adjoint operators at the same time. Now the full invariance group possesses the following Abelian 
subgroup of order 16 

0r = ± {/, Ai, A 2 , As, AiA 2 , AiA 3 , A 2 A 3 , A1A2A3} , (79) 

which is isomorphic to the direct product of cyclic groups of order 2 

g 1 ~ {I, -1} ® {I, Ax} <g> {7, A 2 } <g> {7, A 3 } ~ C 4 (80) 

The elements of Gi will be called dynamical symmetries. With all these Q' can then be expressed as the 
group of products of unitary operators 

g' = {AP\ Aeft.Pefe}. (81) 

This is a finite group of order 16 x 384 = 6144. Development of its representation theory can be based on 
the very important observation that this is in fact a semi-direct product group, because it can be shown 
that Q\ is actually an invariant subgroup. Consequently we can write 

G' = Qi®Q2 ~ C 4 ©0 4 . (82) 

We have already pointed out in the previous section and in particularly Appendix [Cj that this sort of 
group structure is rather convenient when one wants to explore the representations. Convenient, because 
the method of induction |29j can be applied in order to induce irreducible representations of G' from those of 
the constituents Gi and G2, respectively. The method does indeed verify that the irreducible representations 
of Q' have dimensions 1, 4, 6, 12 and 24, so these are to be associated with the energy levels of the confined 4D 
harmonic oscillator (first column of Table [1]). The corresponding noninteracting wave functions of Eq. (|44bj) 
form in turn bases for some (equivalent) forms of these representations. 

4.2. Noninteracting level structure: single- and two-particle density of states 

We now turn to the results of the numerical solution. In Subsection 13.41 we introduced the low energy 
subspace Tl m and constructed a symmetry adapted basis for it. In this every member transforms as some 
row of some irreducible representation. The first few energy levels computed numerically from Eq. (|69[) are 
tabulated in the second column of Table [TJ For numerics we used ui 2 /2 = 500, b = 1, m = 8 and n = 30. 
Consequently the dimension of the whole Hilbert space is n 4 = 810000, whereas that of the subspace is 
to 4 = 4096. The results are represented to six digits of precision and as we can see, compared to the exact 
theoretical values in the first column, remarkable accuracy is found. As a matter of fact this is not that 
surprising as the basis elements of Eq. l|65p are the exact solutions themselves. What is notable though 
is that the computed data are accurate up to three or more digits already at the modest resolution of 
h = 2b/ (n + 1) = 0.0645. This is apparently due to the very small 0(h 6 ) error of our difference scheme. 
This accuracy might give us hope for when we finally incorporate the Coulomb interaction, for which we 
do not have exact analytical results to compare with, the numerical data will be adequate. From the fourth 
column on we can see the distribution of eigenvalues across representations. The last column summarizes 
the degeneracy of a given level. Note that representations written in boldface are two-dimensional. 

So far we have performed the numerical computations at only one step size corresponding to the resolution 
n = 30. In order to verify the conjecture that the eigenvalues do indeed converge fast to the theoretical 
predictions, with a leading error of 0(h e ), we have made test runs with increasing n from n = 10 to 30 with 
a step of 2. Meanwhile, though the dimension of the full Hilbert space grows as n 4 , the dimension of the 
low energy subspace was held fixed as specified above. Figure shows the dependence of the raw data on 
the resolution that were obtained from Eq. (|69| for T 11 . The circles denote the development of the lowest 
energy eigenvalue in this symmetry channel. Similarly, the squares, diamonds and the triangles denote the 
second, fifth and the seventh largest eigenvalues, respectively. As a function of n each of these follows an 
algebraic tendency 

E n = H — r + ... , (83) 
n° 

where stands for the asymptotic limit, i.e. the exact result tabulated in the first column of Table [TJ 
The coefficient a is unknown and also is the exponent b. However, according to what has been said so far 
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n 

Fig. 2. (Color online) Some energy eigenvalues belonging to T 11 as a function of the resolution n. The circles, squares, diamonds 
and triangles are the raw computed data representing the lowest, the second, the fifth and the seventh largest eigenvalue of 
this symmetry, respectively. The solid lines are linear fittings and the legend shows the obtained exponents b, see Eq. J83D in 
the text. Note that both axes are in logarithmic scale. 

we expect b to be close to 6. Linear fittings on the raw data were performed and are shown in the figure 
by solid lines. From these we found that the exponent for the lowest energy (circles) is b — 5.27 ± 0.07, 
for the second largest energy (squares) it is b = 6.29 ± 0.08, for the fifth largest energy (diamonds) it is 
b = 5.95 ± 0.03 and finally for the seventh largest eigenvalue (triangles) b — 5.84 ± 0.04 was obtained. We 
consider these representative results satisfying. The deviations from the theoretically expected exponent 
could be caused for example by the fact that at this relatively low scale of n higher order terms in Eq. (f83|) 
(not shown explicitly) can contribute a little. A more important cause, however, might be the fact that 
we have "truncated" the full Hilbert space considerably, which we did on physical grounds. Increasing the 
dimension of the low energy subspace 7i m would result supposedly in more accurate results and exponents 
closer to 6. 

The level structure of a noninteracting system is usually characterized by the single-particle density of 
states (DOS). In our case, recalling Eq. (|44a[) , it reads 

2 oo 

■9i( e ) -EE^( £ - w Ki + ^ + !)) • ( 84 ) 

i=l m=o 

From this the total noninteracting two-particle density of states can be expressed as a convolution 

4 oo „ 

ff * ot (e) = 6 ( £ - + v »* + "na + "r* + 2)) = / de'gi(e - e') 9l {e'). (85) 

i=l m=o 

This quantity is, however, rather artificial because the spectrum it measures is actually that of an abstract 
four-dimensional particle, the symmetry of its quantum state we do not have restrictions for. Mathematically 
speaking, the Hilbert space of this entity is the direct sum of the symmetric and antisymmetric two-particle 
subspaces. In reality what we have is the pair of two identical particles and the symmetric and antisymmetric 
solutions could describe completely different states. For example consider a spin singlet two-electron state. 
We know that the spatial wave function can only be symmetric in this case, whatever its energy is, and 
so must reside somewhere in the symmetric subspace. As opposed to this the triplet state can only be 
associated with antisymmetric wave functions which then obviously restricts the possible energies of the pair. 
Accordingly, in order to obtain true two-particle DOS, g\ ot has to be split into two terms, loosely speaking, 
with respect to their symmetries under permutation. However, the form in Eq. (|85p is not appropriate for 
this as the quantum numbers rii are not related directly to symmetries. Irreducible representations provide 
good quantum numbers for this purpose and we can write 
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qp r—1 



(86) 



Here E^P is the rth largest eigenvalue of Eq. (|69p . associated with the state ipj p r . The first few noninteracting 
data are shown in the second column of Table [TJ As we have already pointed out, group theory says that 
all energy levels belonging to a certain representation, say T qv , are exactly d 9P -fold degenerate and this 
explains the weight of the dirac-delta. Most importantly, the sum running through representations now 
involves only those for which the basis functions ■0?^. are symmetric (or antisymmetric) under permutation, 
see Subsection 13.51 This constraint is denoted by a prime. 

The insets in the left and right panels of Fig. [JJ show the noninteracting two-particle DOS associated with 
the antisymmetric and symmetric subspaces, respectively. Both gi and g 2 will be modified by the interaction. 
This we will discuss in Subsection 15. II 



4.3. Noninteracting two-particle densities 



Having obtained a particular pair of solutions (E qp , ipj r ) from Eq. (|69p . the corresponding four-dimensional 
wave function (ipj P r ) P iki is then obtained by rearranging the vector index by means of Eqs. (|2T|) and ([22)) . 
From this the two-particle density in the 2D square domain is calculated as 

n 

n pi = J2 ^Ztiki- (87) 

k,l=l 

The lowest energy (r = 1) noninteracting "symmetric" densities belonging to T lp can be seen in the left 
column of Fig. [31 For p = 5 the representation is two-dimensional and the figure shows the density associated 
with the second row, j ' — 2. For j — 1 a very same structure is found that is elongated in direction y instead 
of x. For the other representations we do not show noninteracting densities separately, because they are only 
slightly different from those with interaction depicted in Fig. [U 



5. Numerical results for interacting particles 



In this section we incorporate the effect of a repulsive Coulomb interaction and calculate some relevant 
physical quantities. These are the interacting level structure, one- and two-particle density of states, wave 
functions and particle densities categorized by representations and entanglement. The interaction itself reads 

V{\rx - r 2 |) = C , (88) 

yJ{Xi - XzY + (X 2 - X A y 

where c is a positive constant responsible for its strength. Appearance of this term precludes separation of 
variables in the Schrodinger equation and all dynamical symmetries we encountered during the study of the 
noninteracting problem will be lost. In addition to that many geometrical symmetries will be broken too as 
an arbitrarily chosen R 6 O4 does not necessarily commute with Q of Eq. 146[> . These considerations have 
led to the result that the invariance group Q of the interacting Hamiltonian consists of only geometrical 
symmetries, and as it turned out in Subsection 13.21 it is a subgroup of O4 of order 32. 



5.1. Interacting level structure: single- and two-particle density of states 



Numerical solution of Eq. (|69|) for each irreducible representation and with interaction results in the 
interacting spectrum. For a coupling strength of c = 1 the first twelve levels in each symmetry channel 
are tabulated in Table O In order to put it in a more expressive form, we have also plotted the true 
antisymmetric and symmetric two-particle DOS in the left and right panels of Fig. [SI respectively. For 
evaluation Eq. ([86]) was used with eigenvalues taken from Table [2] The figures reflect clearly that the effect 
of Coulomb interaction in the antisymmetric case is very weak, almost negligible compared to that observed 
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Fig. 3. (Color online) Normalized two-particle densities obtained from the lowest energy symmetric wave functions transforming 
as T lp , p = 1, . . . , 5. In case of p = 5 the row index is j = 2. The left (right) column shows the noninteracting (interacting) 
densities. Associated energy levels in the left (right) column from top to bottom are: 2.000000 (2.702), 4.000019 (4.059), 6.000214 
(6.060), 4.000019 (4.745) and 3.000011 (3.724). These are taken from Tables Q] and [2] 
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Fig. 4. (Color online) Normalized interacting two-particle densities obtained from the lowest energy wave functions. Left (right) 
column shows the antisymmetric (symmetric) densities belonging to T 2p (r 4p ). For two-dimensional representations the row 
index is j = 1. Associated energies in the left column from top to bottom are: 4.078, 4.078, 5.079 and 3.078. In the right column 
from top to bottom: 6.059, 8.048, 6.044, 4.059 and 5.059. These are taken from Tables Q] and [2] 
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Table 2 

Distribution of interacting two-particle energy eigenvalues E'p' /u) across representations of Q. The first twelve levels are shown 
for each representation. The strength of interaction is c = 1 and the other parameters are the same as in Table [T] 
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Fig. 5. (Color online) Two-particle DOS is shown for the antisymmetric (symmetric) subspace in the left (right) panel. Main 
figures show interacting data taken from Table [2] Insets are without interaction. For better visibility a small imaginary part 
(rj = 0.01) was introduced by hand, i.e. the delta-function is approximated by 8(x) « (ti/tt)/(x 2 + r/ 2 ). 

in the symmetric channel. The deviation is easy to understand: the antisymmetry criterion of a state tp P iki 
leads to vanishing of the wave function whenever the two particles are at the same site, VVp» = 0- This in 
turn means that the major contribution from the Coulomb potential is suppressed heavily. In stark contrast 
to this, the diagonal part of the symmetric wave function can be quite large and the interaction, either 
repulsive or attractive, can contribute easily. 

Considering the two-particle problem, g^ is actually nothing else than a nice representation of the level 
structure: positions of the delta-peaks show available bound state energies if we were to introduce two 
interacting particles at the same time in an empty system. The weight of the delta-function on the other 
hand, which is necessarily an integer, accounts for the total number of independent (degenerate) states the 
particles can occupy. Bearing in mind all these, if we go on now with the exploration of the one-particle DOS 
gx, we immediately see that its definition requires a more elaborate analysis. It is because due to interaction, 
introducing a single particle in an empty system or in a system that already has one particle, is not the 
same. Moreover, depending on the symmetry of the two-particle state it can lead to quite different results. 
In the former problem of a single particle living in the system, the interaction obviously does not matter as 
there is no partner to interact with, and g\ is given by Eq. (|84[) . However, if we want to study the spectral 
properties when a second particle is added, we have to apply techniques used in many-body physics. For 
this purpose we introduce the time-ordered zero temperature Green's function |34j 



G(r,r',t) = -i<Xm|T*(r,t)*+(r , ,0)|x, 



(89) 



26 



where \xm) is an occupied one-particle state. Its wave function in coordinate representation and energy 
are, respectively, Xm(r) = ( t>m 1 {x)4>m 2 (y) and £ m = u(v mi + v m2 + 1). Indices rrii can take any nonnegative 
integer values, see Appendix[B] Also, the field operators are most conveniently expressed in this basis: \&(r) = 
Xm(r)a m , where a m is the usual destruction operator. Now this is the point where the discussion has to 
be split into two. If the two-particle state possesses an antisymmetric wave function, like a spin triplet state, 
then a m is fermionic obeying the canonical anticommutation relations. On the other hand if the spatial state 
is symmetric, like in the spin singlet state, then a m has bosonic nature satisfying the usual commutation 
relations. The Lehmann-representation in each case reads 

rW .W^f Wll^^l^)^) , X m (r) Xm (r>) 
( } " \ h h *-E? + e m + iO + e - £m -i0' (90) 

where we have introduced a complete set of orthonormal two-particle states in the first term describing 
particle propagation forward in time. For completeness we note that the second term describes the situation 
when instead of adding an extra particle we remove the particle from its stationary state Xm- Since this is not 
what we want, from now on we will be interested in the first (retarded) term only. Knowing the symmetries 
of the system it is not surprising that we chose the two-particle intermediate states to be stationary states 
of the Schrodinger equation, so they are most conveniently labeled by representations to which they belong. 
Here we would like to call the attention to an interesting property of G, which it does not exhibit in general. 
In true many-body problems the Green's function is usually a very complicated object because neither the 
interacting (ground) state of the TV-particle system, nor the eigenstates of the N ± 1-particle subspaces 
are known. This is why one needs perturbative treatments conventionally. The N — 1 case is, however, 
exceptional. Exceptional because the interacting states are exactly the same as the noninteracting ones. 
Furthermore, the O-particle state is the trivial vacuum state and the 2-particle states are obtained "exactly" 
from numerical computations. Thus everything is given for the construction of G. 

From the knowledge of Green's function g± is expressed in the usual way |34j . In the finite difference 
method, where wave functions are computed only at discrete lattice points, it reads 

•91 W = £' E W r Pd ( 6 ~ E r P + ( 91 ) 
qp r—1 

where the weight of the delta-peak for antisymmetric intermediate states is 

dqp 71 / 71—1 



<-ll' 



= E E E ^ ( Xp )<p n2 (s%)22 2 - ^ (**) ( S XZ ■ ^ 

j— 1 p,i— 1 \ni,n 2 - / 

For symmetric intermediate states it is slightly different 

< = EE E ^i(^)^(^)WC*+^(«p)*m.(«*)WC^ . 03) 

j = lp,i=l \ni,ri2=0 / 

and the matrix element is 

71 

( S Tr)Z"m 2 = E ^nA x p)^n 2 {x i )4)7n 1 { x k)4>77i 2 {xi){^f r ) pl kl- (94) 
p^i,k : l— 1 

The single-particle DOS associated with the antisymmetric and symmetric subspaces are illustrated in the 
left and right panels of Fig. [SJ respectively. The curves belong to that particular case when the first particle 
occupies the ground state with quantum numbers m\ — mi — and energy eo/cu = 2vq + 1 « 1. Insets 
show the respective results without interaction. From these figures the very same conclusion can be drawn 
as from Fig. O Namely, the Coulomb interaction renormalizcs the single particle energies in the symmetric 
subspace more markedly than in the antisymmetric. 
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Fig. 6. (Color online) One-particle DOS is shown for the antisymmetric (symmetric) subspace in the left (right) panel. Quantum 
numbers of the initially placed particle are mi = 7712 = 0. Insets are without interaction. In order to better visualize the 
delta-functions, as in Fig. [5] a small imaginary part (rj = 0.01) was included by hand. 

5.2. Interacting two-particle densities 

The interacting four-dimensional wave functions (4>j r )pikl can be used to compute two-particle densities. 
By means of Eq. (|87p we determined these from the states with lowest energies (r = 1) and plotted them in 
the right columns of Figs. [3] and 2] From Fig. [31 showing densities calculated from symmetric wave functions, 
we can observe again the fact that the Coulomb interaction not only affects the bound state energies but the 
wave functions as well. This in turn leads to strongly modified densities. On the other hand, the interacting 
densities in the left column of Fig. [4l calculated from antisymmetric states, look almost exactly the same as 
without interaction (these are not shown), and this is in accord with our previous finding that the eigenvalues 
were not affected considerably either. 

5.3. Entanglement 

In this subsection we calculate the effect of interaction on the entanglement of the two-particle state. From 
quantum mechanics it is well known that whenever the total wave function tJjab of a bipartite system is 
separable, that is it can be written as a product tpAB — ^a^b-, then the total wave function not only provides 
complete description of the composite system but does it also for both subsystems A and B separately. In this 
case A and B are independent and the terms tpA and ipB are called pure states. In most cases, however, the 
knowledge of the total wave function does not necessarily involve complete characterization of the subsystems 
as a factorization in general does not exist. If this happens, the subsystems are then said to be superposed 
with one another, or in other words are entangled, and each of them is in a mixed state associated with 
reduced density matrices pa and ps, respectively pQ. 

In our case the bipartite system is nothing else than a composite of two identical particles. Each of the two 
is confined in the same 2D square domain so the one-particle Hilbert spaces are the same. Moreover, identity 
of particles imply that though their reduced density matrices themselves may be different, their spectra and 
all other measures of entanglement are the same. Having obtained a particular solution of Eq. (|69[) we can 
construct the density matrix for the "first" particle as 

n 

P P i,p>i' = ^2 W^WW^Wmj (95) 

k,l=l 

where the rows and columns are now labeled by a composite index. In a similar manner, for the "second" 
particle the square of tp is integrated with respect to the coordinates of the first particle. It is well known, and 
can be seen from this expression as well, that p is a nonnegative Hermitian matrix with unit trace jT|. Also, 
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Table 3 

Measures of entanglement calculated from the lowest energy (r = 1) noninteracting and interacting eigenstates belonging to 
five representations of Q. Energy levels in the third row are taken from Tables IT1 and l2l N is the number of Schmidt coefficients 
Aj that are larger than the threshold 10~ e . S is the von Neumann entropy and AS is the difference between values with and 
without interaction. Representations T 11 , T 13 , T 15 and T 42 characterize symmetric states, while T 23 an antisymmetric state. 





r ii 




r 13 




pl5 




p23 




p42 






no int. 


int. 


no int. 


int. 


no int. 


hit . 


no int. 


int. 


no int. 


int. 




2.000000 


2.702 


6.0002f4 


6.060 


3.000011 


3.724 


5.000220 


5.079 


8.002059 


8.048 


Trp 2 


f.000 


0.756 


0.187 


0.140 


0.500 


0.417 


0.250 


0.250 


0.125 


0.082 


N 


1 


30 


6 


8 


2 


34 


4 


8 


8 


11 


S 


0.000 


0.587 


1.733 


2.015 


0.693 


1.079 


1.386 


1.387 


2.079 


2.557 


AS 


0.587 




0.282 




0.386 




0.001 




0.478 





its diagonal elements constitute the usual static particle density we have already encountered in Eq. (|87[) 
and in Figs. [3] and 0J One measure of entanglement is the quantity Trp 2 . This is unity for pure states but 
smaller than one for mixed states, indicating entanglement of the two-particle state. Another characteristic 
measure of superposition comes from the Schmidt decomposition (singular value decomposition) [35j 

N 

Tpviu = Y,^ u< $ v ki- ( 96 ) 

Here u^' and # are, respectively, normalized one-particle eigenstates of the density matrices of the first 
and second particles, each belonging to the same eigenvalue Xj, the so-called Schmidt coefficient. Now, from 
this the von Neumann entropy is found by 

S = - Tip In p = - x j ln A j ■ ■ ( 97 ) 

3 

Both S and the number of nonzero Schmidt coefficients N are widely used measures of entanglement. 

In Table [3] we have tabulated some representative data reflecting the effect of interaction on entanglement. 
In the third row we gave the shift in energy due to interaction for the lowest energy states of five particular 
representations. These data are taken from Tables [1] and [U Further down the table we show the different 
measures of entanglement for noninteracting states as well as for their interacting counterparts. The first 
three and the last representation describe systems with symmetric wave function under permutation of 
particles. The fourth describes an antisymmetric state. As the data verify, the contribution of the repulsive 
Coulomb interaction in the antisymmetric subspace of the two-particle Hilbert space is almost negligible, at 
least for a coupling strength of c = 1. This is in accordance with our previous findings. 

6. Conclusions 

We have studied numerically a Coulomb interacting two-particle system in two-dimensions. This quantum 
mechanical problem called for the solution of the non-relativistic stationary Schrodinger equation, which is 
in fact an eigenvalue equation in a four-dimensional configuration space. As to the domain of motion we 
have specified a square shaped potential well, a computational box with infinitely high walls imposing closed 
boundary condition for the wave function. In addition to that, within the box we have also specified a 
steeply increasing isotropic harmonic potential resembling that in the nucleus. The particles in question are 
identical and can obey either Bose or Fermi statistics. Spin related effects were not considered in this work, 
we concentrated only on the spatial part of the total two-particle state. 

In Section [2] we have developed a fully discretized 89-point finite difference method of the Numerov- 
type that approximates the four-dimensional Laplace operator, and thus the whole Schrodinger equation, 
with a local truncation error of at most 0(h 6 ). The errors of the energy eigenvalues were found to be the 
same order of magnitude. This result, which we have proved analytically as well, can be considered as the 
cornerstone of the present work. It assures that the numerical computations will converge very fast, or, in 
other words, calculations with rather low resolution might give as well reasonable results. We found that the 
finite difference scheme can be put in a remarkably compact and concise form with the usage of matrices with 
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direct polynomial structure. Having obtained the correct matrix representation of the Laplacian we have 
constructed an equivalent generalized matrix Schrodinger equation for the two-particle states and energies. 
In the course of its detailed derivation we came to the conclusion that the required precision, with some 
effort, can be refined further along the same lines. Also, the algorithm can be ported to other dimensions 
too. 

In Section[3]we proceeded with the analysis of internal symmetries of the problem. This we needed because 
the eigenvalue equation, though involves only sparse matrices, is of dimension n 4 and therefore very memory 
consuming. More importantly, we were interested not just in the ground state but also in many excited 
states, and a great deal of unnecessary effort can be eliminated with the use of group representation theory 
and an appropriate similarity transformation that turns the eigenvalue equation into block diagonal form. 
Hence, in this section we determined the invariance group of the interacting Hamiltonian and all of its 
unitary irreducible representations. With the help of these and Wigner-Eckart theorem we performed the 
transformation and obtained a block diagonal eigenvalue equation. Also, irreducible representations allowed 
for a convenient distinction between the symmetric and antisymmetric solutions. 

In Sections S] and [5] we presented the results of our numerical computations. At a resolution of n — 30, i.e. 
with a step size of h w 1/15, the full two-particle Hilbert space is of dimension n 4 = 810000. At this point we 
restricted our studies only to the low energy subspace. This we have chosen to be of dimension 8 4 = 4096. 
In this space we have computed the ground state and roughly the first 200 excited states together with their 
energies. It turned out that in case of the noninteracting results, comparison with exact analytical formulas 
revealed that the numerical data are indeed very accurate, to three or more digits of precision. We classified 
the level structure and the wave functions of the noninteracting as well as the interacting system according 
to the irreducible representations to which they belong. We have computed the static particle densities, 
the one- and two-particle density of states in both the symmetric and antisymmetric subspaces of the full 
Hilbert space. Even investigated briefly the effect of Coulomb interaction on entanglement. To this end we 
have calculated the reduced density matrix, its Schmidt decomposition and the von Neumann entropy. All 
these quantities show consistently that the repulsive interaction affects the antisymmetric states and their 
energies only very mildly, because the vanishing of the wave function when the particles are close to each 
other suppresses the contribution of the Coulomb potential. On the other hand, the symmetric states were 
modified dramatically and the single-particle energies are renormalized considerably as well. 
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Appendix A. Analysis of the matrices Mi, M and N 

In this appendix we summarize some important linear algebraic results on the matrices Mj (i = 1, . . . , 4), 
M and N, that are encountered in Subsection 12.31 There, these properties were used repeatedly during the 
discussion of the discrete Schrodinger equation, Eq. (|3"Tj) . 

A.l. Commutativity 

Taking into account the defining equations l|24p-(|28jl . (|30|) . (|32|) and (|33|) and the properties of Kronecker- 
product, one finds that all these matrices can be put in the common form 



where the indices can take nonnegative integer values. Clearly, in this decomposition it is the coefficients 
Zabcd that are characteristic to the matrices themselves. In case of N they even depend on 7'. The expression 




(A.l) 



a,b,c,d 
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in Eq. (|A. 1|) (written with four different matrices in general) is sometimes referred to as a direct polynomial, 
because it has a multivariable polynomial structure where the variables are matrices and the product is 
the direct-product [19j . Now, knowing the multiplication rule of matrices with direct-product structure, 
remembering that the commutator is a real bilinear function and that every matrix commutes with itself, it 
follows naturally that every two matrices of the set mentioned above commute. In particular 

[M,N] = 0, (A.2) 

and this is independent of 7' and n. A similar result has been obtained in the usual one-dimensional fourth- 
order Numerov method 

A.2. Spectral representation of direct polynomials: the properties o/Mj 

The direct polynomial structure makes it also rather easy to determine the spectrum and the spectral 
representation of such matrices [IS]. With this recognition, next we show that the matrices Mj are all 
negative definite. In our case the basic ingredient is the symmetric n-by-n matrix A, given by Eq. (|27p. 
whose eigenvalues are known analytically 

Wfe = 2cos , ft = 1, ...,n. (A. 3) 

71 + 1 

Since these are all distinct, with algebraic multiplicity of one, A is diagonalizable. Let lit be a normalized 
eigenvector belonging to uik, that is Auj, — uj^k- Then 



1 2 ( . kn 2kn . nkn \ . . 

Tit — \ sin , sm , . . . , sm , (A. 4) 

k Vn+1\ n + 1 ' n + 1 ' n + 1 J ' y ' 

and as is known, the set of these vectors form an orthonormal basis of R™. We note here that a vector u in 
general will always be thought of as a n-by-1 column vector, the respective 1-by-n row form, whenever needed, 
will be explicitly denoted by matrix transpose u T . With these we can write the spectral representation of 
A as 



A = ^2 ^feUfe u fc > ( A - 5 ) 

k=l 

and the corresponding spectral representation of the direct polynomial in Eq. (IA.1|) reads 

H ^abcd^k^k^k^t Ki ® u fc 2 ® u fe3 ® u fe4 ) o (u£ ® u£ $ uj 3 ® uJJ . (A.6) 

Here, k = (fci, k2, ^3, k<±) is a convenient composite index and in the sum each component runs in the range 
1, . . . , n. This is a very useful result for it shows, that though our matrices are rather complex due to their 
multiple direct-product structure, all eigenvalues can be obtained with practically no effort. For example 
the n 4 eigenvalues of Mi are 

mi k ) = 2 cos— — + 2cos^— + 2 cos— 1- 2 cos— 8, h = l, ...,n. (A.7) 

w n + 1 n + 1 n + 1 Ti + l 

With a very same analysis the eigenvalues of M2, M3 and M4 can be given too. They are 

fti7T ki'n k\"K k 3 n k\~K k^K 

wi2 (k) = 4 cos cos h 4 cos cos h 4 cos cos ■ 



and 



71+1 n+1 n+1 n+1 n + 1 n+1 

k'i'K k^TT k2~IT k^TT k^ir k^ir 

+ 4 cos cos h 4 cos cos h 4 cos cos 24, (A. 

n+1 n + 1 n + 1 n + 1 n+1 n+1 

ft!7r k 2 TT k 3 n kxir k 2 ir k 4 ir 
m 3 (k) = 8 cos cos cos + 8 cos cos cos 



71+1 n+1 n+1 n+1 n+1 n + 1 

fcl7T fc37T k^TT klK k^TT k^K 

> cos cos cos h 8 cos cos cos 32, (A. 9) 

n + 1 n + 1 n+1 n + 1 n+1 n+1 v ' 
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and finally 

fcl7T &27T fc 3 7T k^TT 



K\TX K2 71 K 3" K 4 71 

rnAK] = 16 cos cos cos cos 16 

y ' n + 1 n + 1 n + l n+l 



4 cos 2 — — + 4 cos 2 —=— + 4 cos 2 — — + 4 cos 2 — 16. (A.10) 

n+l n+l n+l n+l 



With all these, for any finite n one readily finds 



maxTOi(k)<0, i=l,...,4. (A. 11) 

k 



It shows that the matrices M,; are indeed negative definite. In addition to that 



lim max mi (k) = 0, (A. 12) 

n — >oc k 



that is Mi become singular asymptotically. 



A. 3. The properties of M 



The positivitjQ of M, that we have already adverted to in Subsection 12.31 follows immediately from the 
negativity of Mj and Eq. (|32)) , since all coefficients are negative. Let 8(k) be the eigenvalue of M, then 

0(k) - -i-(12mi(k) + m 2 (k) + m 3 (k)), (A.13) 



30 s 



and its minimum and maximum are 




min k 6>(k)\ _ 16 3 tt 4 2 tt 16 tt 76 



maxk 0(k) 



=F — cos -3 cos 2 =f — cos 1 . (A. 14) 

15 n + 1 5 n+1 5 n + l 15 ; 



Equation (|A.13|) comes from the general result of Eq. (|A.6[) along the same line as mj(k) was obtained in 
Subsection lA.21 The latter result regarding extrema, however, needs some explanation. The allowed k vectors 
are restricted to a cubic domain of M 4 , where their coordinates can only take integer values in the range 
1, . . . , n. In order to find the global minimum and maximum we have to first look for possible local minima 
and maxima inside the domain. The stationary points k* are defined by the solution of Vk#(k*) = 0. Now 
it is easy to show that at these points, wherever they are 

^§^ = 0, i = l,...,4, (A15) 

thus the Hessian matrix composed of the second derivatives is neither positive nor negative definite. This in 
turn leads to the observation that there are no local extrema inside the domain. Consequently, the global 
extrema must be somewhere on the surface of the hypercube. In a similar fashion, again by means of 
Eq. (|A.15|) , one can easily prove that there are no local extrema on lower dimension parts of the surface 
either. Therefore, the global extrema must reside somewhere in the corners. Note that the functions mj(k) 
are symmetric under all permutations of vector components. This feature calls for a classification of the 
16 corners and we find 5 classes with representatives: ki = (1,1,1,1), k 2 = (1,1, l,n), k3 = (1,1, n, n), 
kt = (1, n, n, n) and k5 = (n, n, n, n), respectively. With all these it is sufficient to compare values at class 
representatives: let 9i = 0(kj), then 



A matrix, or more generally a linear operator is called, in brief, positive (negative) if it is positive (negative) definite. 
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16 o 7T 4 , 7T 16 7T 76 , . „ „, 

6i = -— cos J - -cos 2 - — cos + — , (A. 16) 

15 n + 1 5 n+1 5 n + 1 15 v ; 

8 , 7T 8 7T 76 , . „ 

2 = — cos cos 1 , (A. 17) 

15 n + 1 5 n + 1 15 V ; 

6) 3 = Acos 2 ^— + — , (A.18) 

J 15 n + 1 15' _ V ' 

/, 8 o 7T 8 7T 76 , . „ „. 

4 = cos d h-cos 1 , (A. 19) 

15 n+15n+115' y ' 

16 o 7T 4 9 7T 16 7T 76 . . . 

6*5 = — cos cos 2 1 cos 1 . (A. 20) 

5 15 n + 1 5 n+1 5 n+1 15 V ; 

Comparison of these formulas shows that 6> m i n = 6>i and 6* max = 9 5 for all n, which proves Eq. (|A.14|) indeed. 

We end the discussion of M with the presentation of some asymptotic formulas. From Eq. (|A.14|) one 
can easily extract the information that mm and max are strictly decreasing and increasing functions of n, 
respectively. In particular, asymptotic expansion yields 

4^ 2 / 2 



1-- + ... , (A.21) 



"mm — o 

n \ n 

*_-^(.-^ + ...), ,A. 22) 

thus the spectrum obeys 

p(M) = {0 mln , . . . , 6 max } c [0, 128/15]. (A.23) 
The most interesting asymptotic expression, however, is 

( B+ l,. W _ <W ( 1 -|g + £0 +...), (A.24) 

where 

4 

e(k)=^ 2 ^fc 2 , ki = l,2,... (A.25) 

i=l 

is nothing else but the allowed energy levels of a particle in a four-dimensional infinite square well. This 
we have already encountered apropos of Eq. ()34[) in the main text. There we emphasized that these are 
the eigenvalues of the four-dimensional negative Laplace operator supplied by the boundary condition of 
Eq. (|4dp. Equation (|A.24[) thus indicates that (n+ 1) 2 M evolves smoothly into —A, and the error is of order 
(D(n~ 2 ). As we shall see next, this will be refined considerably, to 0(n~ 6 ), by taking into account the matrix 
N as well. 



A. 4. The properties of N and the issue of 7' 

Wrapping up our discussion on matrix properties, we finally turn our attention to the properties of N and 
the determination of 7'. Its eigenvalues A(k) are 

A(k) = 1 + (l2 7 ' - mi(k) + (± - ArA m 2 (k) + 7 'm 3 (k) - ^™ 4 (k). (A.26) 

In Subsection 12.31 apropos of the derivation of the discrete Schrodinger equation and Eq. (|3"4|) , we saw that 
N plays a central role in the 0(h 6 ) theory. In particular, it shows up in the discrete realization of the 
Laplacian. According to arguments given there in 3. and 4., 7' must be chosen so that N is positive definite 
|llj . In other words, even its smallest eigenvalue must be positive. This criterion requires us to calculate the 
global minimum 

A min = minA(k). (A.27) 
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The algorithm we applied in case of M in Subsection IA.3I can again be utilized: stationary points are now 
determined by VkA(k*) = 0, and for the second derivatives we find 

<9 2 A(k*) 7T 2 . 2 k*TT , 

afc/ 30(n + 1)^ Ti + 1 

These are apparently strictly negative, wherever k* may be. The Hessian matrix is surely not positive 
definite, A can have no local minima. The global minimum therefore must be somewhere at the corners of 
the hypercube. First of all, A can be evaluated at the class representatives k^, defined in Subsection IA. 31 Let 
Xi = A(kj), then 

Ai = -h cos4 + 32y cos3 iri + (I - 96y ) cos2 ^ttt + ( 96y - £> cos iri + % - 32y - 

(A.29) 

167' cos 3 — cos 2 — ^— + (487' - — ] cos — ^— + \=- 32 7 ', (A.30) 



A 2 - 


1 4 
COS — 

15 n 


A 3 = 


1 4 
COS 

15 


A 4 = 


1 4 

15 n 


A 5 = 


1 4 
COS 

15 



1 15 n + 1 V 15 / « + 1 15 

J- + f32 7 ' - cos 2 + H _ 327', (A.31) 

n + 1 \ 45/ ti + 1 15 

- + I67' cos 3 - — cos 2 + — - 48 7 ' cos + — - 32 7 ', A.32 

1 ' n+1 15 n+l V 15 J n+1 15 v ; 

7T , o 7T /3 ,\97T / 4 A 11 / 

- - 32 7 ' cos 3 — — + - - 96 7 ' cos 2 — — + — - 96 7 ' cos— - + — - 32 7 ', 

1 7Z. + 1 \o / 71+1 \ 15 / 71 + 1 15 

(A.33) 

and with these we can write 



Amin = minAj. (A. 34) 

i 

Analysis of Xi as functions of 71 and 7' reveals that as long as 

23 

7 ' " 3840 ^ °-°° 59 ' (A - 35) 
Xi are positive for all n. Consequently, the smallest is positive too, resulting finally in the positivity of N. 
Now the ultimate question arises: Is there a specific value in this range whose choice would provide the 
best possible description of the discrete Laplace operator, (71 + 1) 2 N _1 M? The answer is simple: Yes, and 
there is only one! To prove this statement and determine that distinguished value, we have to first study 
the structure of this matrix product. 

Due to commutativity, see Eq. (|A.2|) . it is obvious that they possess common eigenvectors and every 
eigenvalue of N^M is of the form X^ 1 9 : where A and 9 are some eigenvalues of N and M, respectively. 
Moreover, as they exhibit the same direct polynomial structure, we can immediately write 

N^M - TTTT ( u *i ® u k2 ® ufe 3 ® UfcJ « ® u£ ® u£ ® u feJ ■ ( A - 36 ) 
k A[k) 

The spectral representation shows unambiguously that in the product A -1 # mentioned above, both terms 
have to be taken at the very same k. This in turn leads to the fact that every eigenvalue can be obtained in 
such a way, by simply running k over all allowed values. Now, having obtained the whole spectrum, we can 
study its asymptotic behavior for large n. We obtain 

x2 0(k)_^A a(k)-6(k)Y 



where e(k) is given by Eq. (|A.25|) . This is a fundamental result for it shows that the appearance of N _1 in the 
discrete Laplace operator has led to a remarkable refinement of energy levels. Compared to Eq. (|A.24|) we see 
that the order of error dropped to 0(ti -6 ), which means that the numerical procedure converges much faster 
to the exact analytical results. Equation (|A.37[) can be regarded as the analytic proof for that our 0(h e ) 
theory is indeed that accurate. This proof can be considered as a generalization of the results of Refs. [H] 
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and [12) . These papers, based on a matrix norm approach, provide analytical proofs in one-dimension for 
that the local truncation error in the standard Numerov method, which is 0(/i 4 ), determines the accuracy 
of the eigenvalues as well. A more general but from practical point of view less constructive proof is given 
in Ref. [16] for the eigenvalue problem of linear second order self-adjoint and elliptic differential operators 
in arbitrary dimensions. 

Another interesting observation is that the order of error is not affected by 7', there is still a large degree 
of arbitrariness in its specific choice. Though we cannot make it smaller any further, we can still try to 
minimize its prefactor. The right hand side of Eq. (|A.37[) indicates there is, in principle, a perfect choice for 
7' for any given k, the one that would make the leading error vanish. However, detailed calculations show 

k \b(k) J 6(ki) 30240 v ; 

so even the minimum of these values is slightly out of the range specified by Eq. (|A.35|) . From this it follows 
that the most precise approximation we can ever achieve in this theory is the one where 7' is chosen to be 
closest to this minimum, namely 

23 

^ = 3840- (A.39) 
One can also prove the positivity of b(k) for all k, thus Eqs. (|A.37j) and (IA.39j) together imply that the exact 
energy levels are always approached from below. 

Now that we have successfully obtained a very concrete value for 7': can return to the issue of A m i n - 
Comparison of Eqs. (|A.29[) - (|A.33[) reveals that for this choice, A5 is the smallest for all n. Moreover, it can 
also be shown that in this case there are no local maxima of A in the range hi — 1, . . . , n. Hence the global 
maximum must also be among A^, in fact it is A max = Ai. Explicitly 

Amin \ 1 a 7T 23 n 7T 1 9 7T 37 7T 13 , . 

= cos =F cos J 1 cos =F cos 1 . (A. 40) 

\ I 15 n + 1 120 ?i + 1 40 n + 1 120 n + 1 24 ^ ; 

^max / 

From these it turns out that A m i n and A max are strictly decreasing and increasing functions of n, respectively. 
For large arguments 

llTT 2 / 2 \ /A . 

A mi „ = — ^ 1- - + ... , (A.41) 



20 n 2 \ n 
r 2 
3 n 2 



Itt 2 

A nu « = l---a+..., (A.42) 



thus the spectrum obeys 



P (n (3I5)) = {Amin, ■ • ■ , A max } c [0, 1]. (A.43) 

Although N is positive definite, asymptotically it becomes singular at the same rate as M, see Eq. (|A.21|) . 
This is an important observation signalling the fact that the numerical solution of the discrete Schrodinger 
equation in Eq. (|3T|) should not rely on N _1 , because N becomes more and more ill-conditioned as its size 
grows. In fact, its condition number is 

log(^) -logn, (A.44) 

\ A m in / 

which is practically 0(1) in all relevant cases. It is by all means much much larger than the numeric 
(computational) precision of its matrix entries. 

A. 5. Ground state energy of the Laplace operator 

Essentially all expressions and results we obtained in this appendix so far show up in the expression for the 
noninteracting ground state energy. From Eq. (|A.25p we know the exact analytical result: cq = e(ki) = 4-7T . 



35 



Not surprisingly, it is the first energy level of the discrete Laplace operator that is calculated with best 
precision 

•<»> = <" + 1)3 t> (m) = <" + "Hi = <" + 1)2 ^ (A - 45) 

where # m i n and A max are shown in Eqs. (jA.f 61) and (|A.40j) . respectively. Asymptotic expansion yields 

e (n oo) = 4tt 2 fl 703_ttJ A (A 46) 

V ; V 60480 n 6 / V ; 

The prefactor in the error is roughly 10, so even a modest resolution of n — 1 can lead to a relative error 
of 10 -5 , and this is very accurate indeed. 

Appendix B. Harmonic oscillator in a box 

This appendix is devoted to a short study of the one-dimensional quantum harmonic oscillator that is 
spatially confined in a box. We collect some useful results and formulas and also tabulate the first few allowed 
discrete energy levels of the system. Similar studies have been performed in Refs. [26 36 37 3 8139] . 

In Eq. (|43p we defined the dimensionless quadratic potential we use for the numerics throughout the 
paper. In one-dimension 

U(x) = ^u 2 x 2 (B.l) 

for |x| < 6, otherwise U(x) = oo. The parameter to > is a scaling constant responsible for the overall 
strength of the potential. The condition that the oscillator is confined means there is an infinite repulsive 
wall at \x\ > b, which the wave function ip( x ) cannot penetrate into. Within the box the wave equation reads 

-^ + (r v -"H)) ,4 = ' (R2) 

where the energy eigenvalue is written in the form E = uj{v + 1/2). Note that at this stage nothing is known 
about the range of v, certainly except that v > —1/2, as the Hamiltonian is a nonnegative operator. If 
the barrier is removed by the limiting procedure b — > oo, the problem becomes that of the unconstrained 
oscillator. Then the boundary condition is related to the asymptotic behavior of i/j, namely, bound state wave 

functions must be square integrable and this leads finally to the familiar result v = n with n = 0, 1, 2, 

Coming back to Eq. (|B.2[) we see that it is the differential equation of parabolic cylinder functions 
D u , sometimes referred to as Weber equation |40) . Two linearly independent solutions would be natu- 
rally D u {^/ujx) and D v (—^/Zjx). However, these are neither odd nor even functions of x, hence they are not 
suited well for the particular problem. Nevertheless, we can construct appropriate linear combinations that 
are eigenfunctions of parity 

= V^xe-^M Q - V ~, | , (B.3) 

M *) = e-»^M(-l\^), (B.4) 

so ipi t i> is odd and ip2,v is even. Here M(a, 6, z) is the confluent hypergeometric function of the first kind [40) 
defined by the hypergeometric series 

a a(a + 1) z 2 , 
M(a ) M) = l + r+ ^ TI f ¥ + .... (B.5) 

It is an entire function of z provided that 6^0, — 1, —2, .... Now the general solution of Eq. (|B.2[) is 

i> v (x) = A^ hv (x) + BihA*)i ( B - 6 ) 
with A and B being arbitrary constants. 
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Table B.l 

Allowed energy levels of the one-dimensional confined harmonic oscillator, E n = w(u n + 1/2), for the case of b = 1 and 
ui 2 /2 = 500. 



n 




n 







0.000001 


1 


1.000017 


2 


2.000235 


3 


3.001945 


1 


4.010898 


5 


5.043776 


6 


6.132232 


7 


7.315886 


8 


8.628132 


9 


10.088573 


10 


11.705530 


11 


13.481490 


12 


15.416694 


13 


17.510727 


14 


19.763071 


15 


22.173266 



The existence of discrete energy levels is related to the boundary condition imposed by the barrier 

= Mb) = MiAb) + B^ 2 , v [b), (B.7) 
= M-b) = -Aih., v (b) + Bifa t „(b). (B.8) 

Since this is a homogeneous system of linear equations, in order for a nontrivial solution for the coefficients 
A and B to exist, the determinant must vanish. This in turn leads to 

w (1 v 3 u;b 2 \ / v 1 ujb 2 \ n 

Ml ,-, 1 M ( — ,-, =0. B.9 

\2 2' 2' 2 J \ 2 '2 '2 J V ' 

This is the final result that implicitly determines v and such the allowed energy levels. Clearly, this product 
vanishes if either of the two terms becomes zero, thus we have to distinguish two cases: 

(i) Let's consider first 

m(I--4)=o. 

It has infinitely many positive roots which we index with nonnegative odd integers: v\ < v-& < v§ . . . . 
For the particular case of b = 1 and w 2 /2 = 500 the first eight of these are tabulated in Table |B~T1 The 
choice for these parameters is the same we used for numerical computations in Sections [5] and [S] Going 
back to Eq. (|B.7[) or (|B.8|) and setting v = v n it is easy to see that B = must be. The corresponding 
normalized eigenstate is 

M*) = 1n niX n > n= 1,3,5,..., (B.ll) 

where the norm is as usual 

||Vi,,„|| 2 = f iptuM^- (B-12) 

J-b 

It is interesting to examine Eq. (|B.10[) for large v. Asymptotic expansion of M reveals [10] 

W ^ + ^=W (2m + 2)2 ' (R13) 

where m is a large but otherwise arbitrary integer. This shows that at high energy the spectrum 
reproduces half of the familiar energy levels of a particle in a one-dimensional box of size 2b, because 
2m + 2 is always even. 

(ii) Consider now the other equation 

In complete analogy with 1. this has infinitely many positive roots which we now index with nonnegative 
even integers: vq < v 2 < . . . . For the same parameter values as in (i) the first eight of these are 
shown in Table iBTTl Also, from Eq. (IB.7j) we find A — and the normalized eigenstate is 

Mx)= tT ii ' n = 0,2,4,..., (B.15) 
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Table C.l 



Character tables of A (left) and B (right). As A is Abelian all irreducible representations are one dimensional and the characters 
themselves are the matrix elements. Conjugacy classes are denoted by d. 



Ci 


= Ri C 2 = 


R.2 C :i = 


= R 3 C 4 = 


R 4 


Ci 


= Ri C 2 


= {R 4 ,R 8 } C 3 


= R 5 d = 


{R 3 ,R 7 j C b = {R 2 ,R 6 } 




1 


1 


1 


1 


Xb 


1 


1 


1 


1 1 


x'a 


1 


-1 


1 


-1 


x% 


1 


1 


1 


-1 -1 


x A 


1 


-1 


-1 


1 


xl 


1 


-1 


1 


1 -1 


Xa 


1 


1 


-1 


-1 


x% 


1 


-1 


1 


-1 1 












x% 


2 





-2 







where the 


norm 


is analoj 


I to that of Eq. (|B.12|). 


For completeness we 


examine 


Eq. (|B.14| for large v 



too. Asymptotic expansion again provides 

K* + £) = w (2ni+1)2 - (R16) 

This result complements Eq. (|B.13I) as it accounts for that part of the high energy spectrum that is 
related to the square of odd integers. 

Having obtained the allowed energy levels and the corresponding normalized wave functions we can now 
construct the Hilbert space of this one-dimensional problem. As eigenfunctions of the Schrodinger equation 
belonging to different eigenvalues are necessarily orthogonal, the functions <fi n (x) defined by Eqs. (|B.11[) and 
(|B.15[) form an orthonormal set. If the Hilbert space is defined as all possible linear combinations of this set, 
then completeness follows naturally. 

Summing up the results, we found that the infinite potential barrier superimposed on the quadratic 
potential shifts the allowed energy levels upwards, but otherwise, as expected, leaves the discrete feature of 
the spectrum unaffected. At high energy the spectrum and wave functions turn into those of the "particle 
in a box" problem. For the particular choice of parameters we found that close to the ground state the wave 
functions and their energies are practically those of the unconstrained oscillator. Nevertheless, in higher 
dimensions these shifts, whatever small they are, will lift certain degeneracies of the excited energy levels 
and so lead to qualitative changes. 

Appendix C. Irreducible representations of Q 

In this appendix we explicitly construct all inequivalent unitary irreducible representations of the group 
of the Schrodinger equation. In Eq. (|50|) we have already pointed out that Q has a semi-direct product 
structure. This is a rather satisfactory situation, because in this special case the knowledge of all irreducible 
representations of A and B is sufficient to induce those of Q by means of general theorems of group theory. 
The method itself is called induction. Before that, however, we shall give a very brief overview of the 
representations of A and B. Representations and characters will be denoted by T and %, respectively. 

C.l. Representations of A and B 

Using the matrices of Eq. |47|) direct calculations may verify that A is indeed Abelian. Irreducible rep- 
resentations are necessarily one-dimensional. Furthermore, as R2 = R3R4 and the square of any member 
equals identity, it has a direct-product structure 

^-{R 1 ,R 3 }®{Ri,R 4 }-C 2 2 . (C.l) 

It shows that A is isomorphic to Cf , with C2 being the cyclic group of order 2. From this the character 
table follows naturally, see Table IC.ll 

The matrix group B was defined by its members in Eq. (jlS)) . As noted there, it is essentially nothing 
else than a faithful four-dimensional representation of the point symmetry group of a square, C^ v . The fact 
it is a group of order 8 and has 5 classes involves it has 5 irreducible representations that are unique up 
to relabelling. Its character system is shown in Table IC.ll The first four are one-dimensional, completely 



38 



described by their characters. The fifth is two-dimensional and a concrete realization with unitary matrices 






r^(R 7 ) = - 



-l 

l 




(C.2) 

Comparison with Eq. (|4"5| shows that these are nothing else but the upper left (or bottom right) 2-by-2 
blocks of Ri, so r| is a faithful representation. 

C.2. Induced representations of ' Q = A(§)B 

In order to obtain the representations r of Q we have to first explore the little groups of B. In what follows 
we apply the same notation as that of Ref. [29] . Let B(q) be the subset of elements R& of B such that 

x^RbRaR,- 1 ) = x q A CR a ) (C.3) 

for all R a £ A. Then B(q) is a subgroup of B and it is called the qth little group. Straightforward calculation 
yields 

B{l) = B(4) = B, (C.4) 
B(2)=B(3) = {R 1 ,R 4 ,R 5 ,R 8 }. (C.5) 

For explicit forms of matrices see Eq. (gHl). According to Eq. (|C.4j) the little groups of q = 1 and 4 are B 
itself thus the orbit of q = 1 is the set {1}, whereas that of q = 4 is {4}. Irreducible representations of Q 
belonging to these values are found easily 29J 

r^(R a R fe ) = X q A (R a )T p B (K b ), g-1,4, p=l,...,5. (C.6) 

Essentially only p = 5 is a real matrix representation because it is two-dimensional, see Eq. (|C.2[) . The others 
are one-dimensional and their matrix elements are tabulated in Table IC.ll 

On the other hand it turns out that the orbit of q — 2 as well as q = 3 is {2,3}. This is the case where 
the method of induction is actually made use of. As it is sufficient to consider only one element in each 
orbit, we choose q = 2. We will need shortly the coset representatives for the decomposition of B into right 
cosets with respect to B(2): they are Ri and R7 of B, respectively. Also, inspection of 23(2) reveals that it 
is Abelian and isomorphic to A 

Ri 6 B < — ► R x e A, (C.7) 
R4 G B < — ► R 4 G A, (C.8) 
R 5 G B < — ► R 3 G A, (C.9) 
R 8 e8< — ► R 2 e A (CIO) 

Because of this the representations of £>(2) are the same as those of A shown in Table IC.ll With all these 
information at hand we are now in a position to explicitly give the matrix elements of the remaining four 
two-dimensional irreducible representations of Q. The upper left elements read 

ifR b GS(2), =i 4 

ifR b £B(2), P 



r 2 p(R a R fc )n = |^ (Ra)xPB(2)(Rb 



whereas the upper right are 



1 ' |0, if R5R7 1 i B(2), p v ' 
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The lower left are given by 



r^(R Q R,) 21 J^(W B(2) (R^ ifR 7 R b e B (2), 

V ' 10, ifR 7 R fc ^(2), * k ' 



and finally the lower right are obtained as 



p 2 p fR R s = fx^(Ra)xg (2) (R 7 R6R7 1 ), if RrRbRf 1 e 6(2), 
° [0, ifRyRbRf 1 £B{2), 



p = l,...,4. (C.14) 



Perhaps the most important irreducible representation of all, though mathematically trivial, is the so- 
called completely symmetric representation given by 

r n (R) = l, for all Re^. (C.15) 

Summarizing the results achieved in this appendix we can say there are altogether 5 + 5 + 4 = 14 
inequivalent unitary irreducible representations of the group of the Schrodinger equation. Because of their 
intrinsic product structure they are conveniently labeled by a composite index qp (no multiplication). If 
q = 1 or 4 then p = 1, . . . , 5, while if q = 2 then p = 1, . . . , 4. Among the fourteen representations there 
are six two-dimensional: qp = 15, 45, 21, 22, 23 and 24, respectively. All others are one-dimensional. Since 
all representations involve only real matrices the unitary property is equivalent to orthogonality. Further, 
from the character tables and the formulas above it is apparent that any matrix element can only be 0, 
1 or —1. This, in conjunction with orthogonality results in the fact that each T 9P (R) is a so-called signed 
permutation matrices: there is exactly one nonzero entry in each row and column and these are either 1 or 
-1. 



Appendix D. Symmetry adapted basis of the Hilbert space 



In Subsection 13.41 we introduced a low-energy subspace TL m in the full Hilbert space. This is spanned by 
all real linear combinations of the vectors v(k) given by Eq. (f65|) . In this appendix we show that the method 
of projections results in a new orthogonal basis in it, every member of which transforms as some row of 
some irreducible representation of Q . This newly formed set is called the symmetry adapted basis and will 
be used in the solution of the eigenvalue problem. 

The four components of k take integer values in the range 0, . . . , m — 1 < n, there are thus altogether m 4 
such vectors. Let us denote the set of all k by K. As to the scalar product we find 

™ 4 

(v(k),v(k')) = 5> M (k)7v(k') = <W, (D.l) 

indicating these vectors are indeed orthonormal. 

Consider next the group of the Schrodinger equation. From either Eqs. (|47|) . (|48|) and (|49|) , or from the 
fact that Q is a subgroup of O4 it follows 

R = d(R)<x(R) (D.2) 

for all R G Q . Here d is a diagonal matrix with ±1 in the diagonal and er is a permutation matrix with 
only one nonzero element in each row and column, which is 1. Actually, <x is nothing else than |R| where 
the absolute value should be taken element-wise. This decomposition is unique and the order of terms is 
important. Writing the terms in opposite order yields 

erd = diag (erdiag(d)) <x 7^ der, (D.3) 

where diag(a) is either a diagonal matrix composed of the column vector a or a column vector extracted from 
the diagonal matrix a. With the aid of this equation it is now not too difficult to verify that the mapping 
<t(R) from Q is a four-to-one homomorphism onto a group of permutations Q* of order 8. Applying the 
scalar transformation operator of Eq. (|56|) to a basis vector we get 

P(R)v(k) = /(R, k)v(o-(R)k), (D.4) 
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where the coefficient is 

4/4 \ fe - 

/(R,k)=n e^O ■ (D - 5) 
1=1 \p=i j 

Careful inspection justifies the conjecture that / is also a homomorphic mapping of Q 

/(R'R, k) = J(R', o-(R)k)/(R, k). (D.6) 

This time it is sixteen-to-one and the mapping is onto the group {1,-1}. Going back to Eq. (|D.4[) we see 
that the basis vectors have a characteristic feature: under coordinate transformations, apart from potential 
sign changes, they transform among themselves. This indicates that Tl m is an invariant subspace for all 
P(R). 

Next we define certain subsets of K. Let K (k) be the subset including the following permutations 

K(k) — {erk | er e Q*} , (D.7) 

and an element is contained only once. It is obvious that every element of K is in exactly one such class 
and the set of all different classes forms a complete disjoint decomposition of K. Combinatorial calculation 
gives for the total number of different classes 

/ m \ I m \ I Tn \ 

(D.8) 

Let us denote the ith class by Ki and its order by Zi. Thus we obviously have 

N 

|J Ki = K, (D.9) 

and 

N 

^z^m 4 . (D.10) 

4=1 

Now apply the projection operators as prescribed in Eq. ([SS]) . Using Eq. (|D.4[1 we get 

wf (k) = Pgv(k) = (d qp /g) ]T r^(R) ii /(R,k)v( CT (R)k), (D.ll) 

Res 

showing that for any given k the new function is either zero or a nonzero linear combination of the original 
basis elements belonging to the class if (k). This means that not just the whole space associated with 
K but also each smaller dimensional subspace associated with Ki are invariant subspaces of the scalar 
transformation operators as well as the projections. For the moment let us assume that the vector in 
Eq. (|D.11[) docs not vanish. Therefore, for any other k' that is not contained in K(k) the resulting wj p (k'), 
if not zero, must be necessarily orthogonal to wj p (k) and so linearly independent. On the other hand, should 
it be any other member of the same class we would obtain 

wf (k') = Pg v(k') = /(R*, k) £ r*»(R*)*P£v(k), (D.12) 

j=i 

where we used that /, er and the representation T qp are all homomorphisms and R* is any element of Q such 
that k' = er(R*)k. Now it is easy to see that if the representation is one-dimensional, that is d qv = i = j = 1, 
this vector is proportional to wj p (k), which in turn means linear dependence. 

If d qp is greater than one the analysis of independence is not so trivial and we can only quote here 
the result: there might be more than one linearly independent vectors, but in any case they are mutually 
orthogonal. Though the proof of orthogonality needs some extra effort, the number of these functions can be 
obtained quite easily as follows. Introduce the quantity rf , which is by definition an integer and measures 
the number of linearly independent functions transforming as, say the jth row of T qp and result from the 
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projections applied on the subspace associated with JQ. They are wj p (k s ), s = 1, . . . , rf. Let us now rewrite 
Eq. CEO) as 

P(R)v(k) = E r(RV k v(k'), (D.13) 

where k S Ki as well and 

r(R) k , k = /(R,k)5 k ,, er(R)k . (D.14) 

As the notation suggests, this is a representation of Q of dimension Zi. Also, because of Eq. (jD.ip the repre- 
sentation consists of orthogonal matrices only. This, in conjunction with Eq. (|D.14[) shows that the matrices 
are actually signed permutation matrices. The number of times rf that an irreducible representation T qp 
appears in it is given explicitly by the respective characters 



rf = - E X,(R)X 9P (R), (D.15) 
5 ReO 

where x qp 1S the character of T qp (see Appendix ICj) and Xi(R-) — J2keK F(R) kk . These also obey 



IP' i 



(D.16) 



A very interesting result is found finally: in any linearly independent set |w' p (k)} (here the indices are fixed 
and k varies) the members are mutually orthogonal. Remember that independence of vectors in inner product 
spaces is an obvious consequence of orthogonality, but the reverse statement does not hold usually. Now 
that we have explored these characteristics we can give an efficient algorithm for the explicit construction 
of the symmetry adapted basis. This method is also easy to implement numerically. 

(i) Consider a class Ki of K and a projection Pj? with indices fixed. Apply it to the elements of the 
class one by one. According to the results above, in the newly formed set {w* p (k)}, which is at most 
of order Zi, there is at most one linearly independent vector if the representation is one-dimensional. 
In this case one should keep the first and forget about the rest. Note, that it might happen that the 
projection results in zero vectors only. In this case go to (ii) because we don't need them: members of 
a basis can only be nonzero vectors. On the other hand, if a representation of dimension greater than 
one is considered the results above tell that there are exactly rf independent functions in the set and 
these are orthogonal. As a result the selection mechanism of independence can be equally based on 
orthogonality which is much faster numerically, at least when there are many (and large) vectors to 
deal with. 

(ii) Repeat (i) for every N classes of K, all irreducible representations qp and all rows j. 

The vectors kept in this algorithm necessarily form a maximal linearly independent set, all of which are 
members of 7Y m . The crucial observation is now the fact that their number is exactly m 4 , which follows from 
the completeness of projections [33] 
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qp j'=i 



Pg = E B *. (D.17) 



This means that the symmetry adapted vectors form orthogonal but not necessarily normalized basis for TL m 
and this is what we wanted to prove. With the aid of Eq. (|D.15|) the total number of independent elements 
in the new basis transforming as say the jth row of T qp reads 



JV 

J-rf. 

1=1 



(D.18) 



The vectors themselves will be denoted by w| p (k s ), s = 1,.. . ,r qp . Similarly to Eq. (|D.16|) the dimensions 
satisfy an analogue constraint 

qp 
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which follows immediately from Eq. (|D.10[) and expresses the fact that the whole subspace associated with 
K decomposes into a direct sum of smaller dimensional subspaces, each belonging to a specific irreducible 
representation. 
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